Evolutionary rescue experiment: effect of genetic diversity
Laure Olazcuaga
07-March-2022
1 DATA
1.1 Pop sive over time dataset
data <- read.table(file=here::here("data", "BE_Census_Data.csv"), sep=",", header=TRUE)
head(data)## Block Old_Label Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High1 High1 High
## 2 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High2 High2 High
## 3 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High3 High3 High
## 4 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High4 High4 High
## 5 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High5 High5 High
## 6 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High6 High6 High
## Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1 0 1 2/17/21 2/17/21 2/18/21
## 2 0 1 2/17/21 2/17/21 2/18/21
## 3 0 1 2/17/21 2/17/21 2/18/21
## 4 0 1 2/17/21 2/17/21 2/18/21
## 5 0 1 2/17/21 2/17/21 2/18/21
## 6 0 1 2/17/21 2/17/21 2/18/21
## Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1 DSC_0155 <NA> NA NA 101.9 NA NA
## 2 DSC_0156 <NA> NA NA 104.7 NA NA
## 3 DSC_0157 <NA> NA NA 105.8 NA NA
## 4 DSC_0158 <NA> NA NA 101.9 NA NA
## 5 DSC_0159 <NA> NA NA 102.2 NA NA
## 6 DSC_0160 <NA> NA NA 101.7 NA NA
## HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1 100 NA NA no no
## 2 100 NA NA no no
## 3 100 NA NA no no
## 4 100 NA NA no no
## 5 100 NA NA no no
## 6 100 NA NA no no
## Less_than_5
## 1 no
## 2 no
## 3 no
## 4 no
## 5 no
## 6 notail(data)## Block Old_Label Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19 Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20 Med20
## 573 Block5 B5-Backup Fam A 8/25 BE B5 | G6 Med 21 Med21
## 574 Block5 B5-Backup Fam B 8/25 BE B5 | G6 Med 22 Med22
## 575 Block5 B5-Backup Fam C 8/25 BE B5 | G6 Med 23 Med23
## 576 Block5 B5-Backup Fam D 8/25 BE B5 | G6 Med 24 Med24
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571 Med 5 6 8/25/21
## 572 Med 5 6 8/25/21
## 573 Med 5 6 8/25/21
## 574 Med 5 6 8/25/21
## 575 Med 5 6 8/25/21
## 576 Med 5 6 8/25/21
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2
## 571 8/25/21 8/26/21 DSC_0273 DSC_0274 2.00000 2.00000
## 572 8/25/21 8/26/21 <NA> <NA> NA NA
## 573 8/25/21 8/26/21 DSC_0277 DSC_0278 29.76978 32.92378
## 574 8/25/21 8/26/21 <NA> <NA> NA NA
## 575 8/25/21 8/26/21 DSC_0267 DSC_0268 25.31021 46.10092
## 576 8/25/21 8/26/21 DSC_0283 DSC_0284 71.09684 55.37443
## MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571 4.0 2 2 4 8 0.500000
## 572 0.0 NA NA NA 0 NA
## 573 62.7 23 23 46 35 1.314286
## 574 0.0 NA NA NA NA NA
## 575 71.4 20 36 56 40 1.400000
## 576 126.5 76 57 133 41 3.243902
## First_Throw Extinction_finalstatus Less_than_5
## 571 no no yes
## 572 extinct yes yes
## 573 no no no
## 574 extinct yes yes
## 575 no no no
## 576 no no nodim(data)## [1] 576 23summary(data)## Block Old_Label Label Population
## Length:576 Length:576 Length:576 Length:576
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## Length:576 Min. :0.0 Min. :1.0 Length:576
## Class :character 1st Qu.:1.0 1st Qu.:2.0 Class :character
## Mode :character Median :2.5 Median :3.5 Mode :character
## Mean :2.5 Mean :3.5
## 3rd Qu.:4.0 3rd Qu.:5.0
## Max. :5.0 Max. :6.0
##
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2
## Length:576 Length:576 Length:576 Length:576
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
##
## MC_Box1 MC_Box2 MC_Total_Adults HC_Box1
## Min. : 0.00 Min. : 0.00 Min. : 0.00 Min. : 0.00
## 1st Qu.: 26.52 1st Qu.: 25.99 1st Qu.: 51.05 1st Qu.: 22.00
## Median : 58.87 Median : 54.29 Median :101.50 Median : 56.50
## Mean : 75.88 Mean : 70.57 Mean :128.37 Mean : 72.06
## 3rd Qu.:106.93 3rd Qu.: 99.40 3rd Qu.:170.25 3rd Qu.:101.00
## Max. :327.40 Max. :299.00 Max. :514.40 Max. :288.00
## NA's :142 NA's :141 NA's :5 NA's :144
## HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## Min. : 0.00 Min. : 0.0 Min. : 0.00 Min. :0.0000
## 1st Qu.: 25.00 1st Qu.: 62.0 1st Qu.: 67.25 1st Qu.:0.4183
## Median : 52.00 Median :100.0 Median :100.00 Median :0.9386
## Mean : 68.08 Mean :132.1 Mean :138.30 Mean :1.4482
## 3rd Qu.: 96.75 3rd Qu.:169.0 3rd Qu.:188.00 3rd Qu.:2.3587
## Max. :290.00 Max. :499.0 Max. :499.00 Max. :6.8571
## NA's :142 NA's :44 NA's :122 NA's :142
## First_Throw Extinction_finalstatus Less_than_5
## Length:576 Length:576 Length:576
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
## #Remove populations killed due to the pathogen
data_status <- data.frame(data$Population,data$Extinction_finalstatus)
data_kill<-data[data$Extinction_finalstatus=="kill",]
dim(data_kill)## [1] 72 23data<-data[data$Extinction_finalstatus!="kill",]
#Add pop growth rate
data$Nb_adults <- data$HC_Total_Adults
data$Lambda <- data$Nb_adults / data$Nb_parents
data$obs<-as.factor(1:nrow(data))
#summary(data)
#Remove
#Check variables
str(data)## 'data.frame': 504 obs. of 26 variables:
## $ Block : chr "Block3" "Block3" "Block3" "Block3" ...
## $ Old_Label : chr "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
## $ Label : chr "BE B3 2/17 | G1 High1" "BE B3 2/17 | G1 High2" "BE B3 2/17 | G1 High3" "BE B3 2/17 | G1 High4" ...
## $ Population : chr "High1" "High2" "High3" "High4" ...
## $ Genetic_Diversity : chr "High" "High" "High" "High" ...
## $ Generation_Parents : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Generation_Eggs : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Date_Census : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_Start_Eggs : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_End_Eggs : chr "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
## $ Image_Box1 : chr "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
## $ Image_Box2 : chr NA NA NA NA ...
## $ MC_Box1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Box2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Total_Adults : num 102 105 106 102 102 ...
## $ HC_Box1 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Box2 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Total_Adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Nb_parents : num NA NA NA NA NA NA NA NA NA NA ...
## $ Growth_rate : num NA NA NA NA NA NA NA NA NA NA ...
## $ First_Throw : chr "no" "no" "no" "no" ...
## $ Extinction_finalstatus: chr "no" "no" "no" "no" ...
## $ Less_than_5 : chr "no" "no" "no" "no" ...
## $ Nb_adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Lambda : num NA NA NA NA NA NA NA NA NA NA ...
## $ obs : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...data$Genetic_Diversity <- as.factor(data$Genetic_Diversity)
data$Block <- as.factor(data$Block)
data$Population <- as.factor(data$Population)
data$Extinction_finalstatus <- as.factor(data$Extinction_finalstatus)
#Order levels
data$Genetic_Diversity <- plyr::revalue(data$Genetic_Diversity, c("Med"="Medium"))
data$Genetic_Diversity <- factor(data$Genetic_Diversity, levels = c("High", "Medium", "Low"))
levels(data$Genetic_Diversity)## [1] "High" "Medium" "Low"data<-droplevels(data) #remove previous levels
str(data)## 'data.frame': 504 obs. of 26 variables:
## $ Block : Factor w/ 3 levels "Block3","Block4",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ Old_Label : chr "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " "B3_All_Red_Mix " ...
## $ Label : chr "BE B3 2/17 | G1 High1" "BE B3 2/17 | G1 High2" "BE B3 2/17 | G1 High3" "BE B3 2/17 | G1 High4" ...
## $ Population : Factor w/ 84 levels "High1","High13",..: 1 9 19 24 25 26 27 28 39 49 ...
## $ Genetic_Diversity : Factor w/ 3 levels "High","Medium",..: 1 1 1 1 1 1 1 3 3 3 ...
## $ Generation_Parents : int 0 0 0 0 0 0 0 0 0 0 ...
## $ Generation_Eggs : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Date_Census : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_Start_Eggs : chr "2/17/21" "2/17/21" "2/17/21" "2/17/21" ...
## $ Date_End_Eggs : chr "2/18/21" "2/18/21" "2/18/21" "2/18/21" ...
## $ Image_Box1 : chr "DSC_0155" "DSC_0156" "DSC_0157" "DSC_0158" ...
## $ Image_Box2 : chr NA NA NA NA ...
## $ MC_Box1 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Box2 : num NA NA NA NA NA NA NA NA NA NA ...
## $ MC_Total_Adults : num 102 105 106 102 102 ...
## $ HC_Box1 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Box2 : int NA NA NA NA NA NA NA NA NA NA ...
## $ HC_Total_Adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Nb_parents : num NA NA NA NA NA NA NA NA NA NA ...
## $ Growth_rate : num NA NA NA NA NA NA NA NA NA NA ...
## $ First_Throw : chr "no" "no" "no" "no" ...
## $ Extinction_finalstatus: Factor w/ 2 levels "no","yes": 1 1 1 1 1 1 1 1 1 1 ...
## $ Less_than_5 : chr "no" "no" "no" "no" ...
## $ Nb_adults : int 100 100 100 100 100 100 100 100 100 100 ...
## $ Lambda : num NA NA NA NA NA NA NA NA NA NA ...
## $ obs : Factor w/ 504 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...head(data)## Block Old_Label Label Population Genetic_Diversity
## 1 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High1 High1 High
## 2 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High2 High2 High
## 3 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High3 High3 High
## 4 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High4 High4 High
## 5 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High5 High5 High
## 6 Block3 B3_All_Red_Mix BE B3 2/17 | G1 High6 High6 High
## Generation_Parents Generation_Eggs Date_Census Date_Start_Eggs Date_End_Eggs
## 1 0 1 2/17/21 2/17/21 2/18/21
## 2 0 1 2/17/21 2/17/21 2/18/21
## 3 0 1 2/17/21 2/17/21 2/18/21
## 4 0 1 2/17/21 2/17/21 2/18/21
## 5 0 1 2/17/21 2/17/21 2/18/21
## 6 0 1 2/17/21 2/17/21 2/18/21
## Image_Box1 Image_Box2 MC_Box1 MC_Box2 MC_Total_Adults HC_Box1 HC_Box2
## 1 DSC_0155 <NA> NA NA 101.9 NA NA
## 2 DSC_0156 <NA> NA NA 104.7 NA NA
## 3 DSC_0157 <NA> NA NA 105.8 NA NA
## 4 DSC_0158 <NA> NA NA 101.9 NA NA
## 5 DSC_0159 <NA> NA NA 102.2 NA NA
## 6 DSC_0160 <NA> NA NA 101.7 NA NA
## HC_Total_Adults Nb_parents Growth_rate First_Throw Extinction_finalstatus
## 1 100 NA NA no no
## 2 100 NA NA no no
## 3 100 NA NA no no
## 4 100 NA NA no no
## 5 100 NA NA no no
## 6 100 NA NA no no
## Less_than_5 Nb_adults Lambda obs
## 1 no 100 NA 1
## 2 no 100 NA 2
## 3 no 100 NA 3
## 4 no 100 NA 4
## 5 no 100 NA 5
## 6 no 100 NA 6tail(data)## Block Old_Label Label Population
## 571 Block5 B5 - Backup Fam H 8/25 BE B5 | G6 Med 19 Med19
## 572 Block5 B5 - Backup Fam I 8/25 BE B5 | G6 Med 20 Med20
## 573 Block5 B5-Backup Fam A 8/25 BE B5 | G6 Med 21 Med21
## 574 Block5 B5-Backup Fam B 8/25 BE B5 | G6 Med 22 Med22
## 575 Block5 B5-Backup Fam C 8/25 BE B5 | G6 Med 23 Med23
## 576 Block5 B5-Backup Fam D 8/25 BE B5 | G6 Med 24 Med24
## Genetic_Diversity Generation_Parents Generation_Eggs Date_Census
## 571 Medium 5 6 8/25/21
## 572 Medium 5 6 8/25/21
## 573 Medium 5 6 8/25/21
## 574 Medium 5 6 8/25/21
## 575 Medium 5 6 8/25/21
## 576 Medium 5 6 8/25/21
## Date_Start_Eggs Date_End_Eggs Image_Box1 Image_Box2 MC_Box1 MC_Box2
## 571 8/25/21 8/26/21 DSC_0273 DSC_0274 2.00000 2.00000
## 572 8/25/21 8/26/21 <NA> <NA> NA NA
## 573 8/25/21 8/26/21 DSC_0277 DSC_0278 29.76978 32.92378
## 574 8/25/21 8/26/21 <NA> <NA> NA NA
## 575 8/25/21 8/26/21 DSC_0267 DSC_0268 25.31021 46.10092
## 576 8/25/21 8/26/21 DSC_0283 DSC_0284 71.09684 55.37443
## MC_Total_Adults HC_Box1 HC_Box2 HC_Total_Adults Nb_parents Growth_rate
## 571 4.0 2 2 4 8 0.500000
## 572 0.0 NA NA NA 0 NA
## 573 62.7 23 23 46 35 1.314286
## 574 0.0 NA NA NA NA NA
## 575 71.4 20 36 56 40 1.400000
## 576 126.5 76 57 133 41 3.243902
## First_Throw Extinction_finalstatus Less_than_5 Nb_adults Lambda obs
## 571 no no yes 4 0.500000 499
## 572 extinct yes yes NA NA 500
## 573 no no no 46 1.314286 501
## 574 extinct yes yes NA NA 502
## 575 no no no 56 1.400000 503
## 576 no no no 133 3.243902 5041.2 Remaining heterozygosity
#Upload growth rate end of experiment
data_he <- read.table(file=here::here("data", "He_ReMaining_BeatricePop.csv"), sep=",", header=TRUE)
data_he$Population <- as.factor(data_he$Population)
data_he <- data_he[,c(1,3)]
### Additional: He remaining
# New vector for the lost during exp
data_he$He_remain_exp <- NA
# He lost during exp
for(i in levels(data$Population)){
temp_data_pop <- subset(data, Population==i)
ifelse(sum(temp_data_pop$Generation_Parents)!=15,print("ERROR"),print(i))
temp_he_remaining_gen <- 1 - (1/(2*temp_data_pop$Nb_parents))
#To consider extinct add this row:
is.na(temp_he_remaining_gen) <- sapply (temp_he_remaining_gen, is.infinite)
temp_he_remain_exp <- prod(temp_he_remaining_gen, na.rm = TRUE)
data_he$He_remain_exp[data_he$Population==i] <- temp_he_remain_exp
rm(temp_he_remain_exp,temp_he_remaining_gen,temp_data_pop)
}## [1] "High1"
## [1] "High13"
## [1] "High14"
## [1] "High15"
## [1] "High16"
## [1] "High17"
## [1] "High18"
## [1] "High19"
## [1] "High2"
## [1] "High20"
## [1] "High21"
## [1] "High22"
## [1] "High23"
## [1] "High24"
## [1] "High25"
## [1] "High26"
## [1] "High27"
## [1] "High28"
## [1] "High3"
## [1] "High32"
## [1] "High33"
## [1] "High34"
## [1] "High35"
## [1] "High4"
## [1] "High5"
## [1] "High6"
## [1] "High7"
## [1] "Low1"
## [1] "Low10"
## [1] "Low11"
## [1] "Low12"
## [1] "Low13"
## [1] "Low14"
## [1] "Low15"
## [1] "Low16"
## [1] "Low17"
## [1] "Low18"
## [1] "Low19"
## [1] "Low2"
## [1] "Low20"
## [1] "Low21"
## [1] "Low22"
## [1] "Low24"
## [1] "Low25"
## [1] "Low26"
## [1] "Low27"
## [1] "Low28"
## [1] "Low29"
## [1] "Low3"
## [1] "Low30"
## [1] "Low31"
## [1] "Low32"
## [1] "Low33"
## [1] "Low34"
## [1] "Low35"
## [1] "Low36"
## [1] "Low4"
## [1] "Low5"
## [1] "Low7"
## [1] "Low8"
## [1] "Low9"
## [1] "Med1"
## [1] "Med10"
## [1] "Med11"
## [1] "Med12"
## [1] "Med13"
## [1] "Med14"
## [1] "Med15"
## [1] "Med17"
## [1] "Med18"
## [1] "Med19"
## [1] "Med2"
## [1] "Med20"
## [1] "Med21"
## [1] "Med22"
## [1] "Med23"
## [1] "Med24"
## [1] "Med3"
## [1] "Med4"
## [1] "Med5"
## [1] "Med6"
## [1] "Med7"
## [1] "Med8"
## [1] "Med9"# Remaining He at the end of the experiment
data_he$He_remain_end <- data_he$He_remain * data_he$He_remain_exp
# Remove kill population
data_he <- data_he[!is.na(data_he$He_remain_exp),]
data_he <- droplevels(data_he)1.3 Phenotyping dataset
#Upload growth rate end of experiment
data_phenotyping <- read.table(file=here::here("data", "Adaptation_Data.csv"), sep=",", header=TRUE)
#Factors variables
data_phenotyping$ID_Rep <- as.factor(data_phenotyping$ID_Rep)
data_phenotyping$Week <- as.factor(data_phenotyping$Week)
data_phenotyping$Block <- as.factor(data_phenotyping$Block)
data_phenotyping$Population <- as.factor(data_phenotyping$Population)
#Revalue
data_phenotyping$Genetic_Diversity <- as.factor(data_phenotyping$Divsersity)
data_phenotyping$Genetic_Diversity <- plyr::revalue(data_phenotyping$Genetic_Diversity, c("Med"="Medium"))
data_phenotyping$Genetic_Diversity <- factor(data_phenotyping$Genetic_Diversity, levels = c("High", "Medium", "Low"))
#Add sum total
data_phenotyping$Nb_ttx <- rowSums(data_phenotyping[,11:13], na.rm=TRUE)
#Proportion
data_phenotyping$p_adults <- data_phenotyping$Adults / data_phenotyping$Nb_ttx
data_phenotyping$p_pupae <- data_phenotyping$Pupae / data_phenotyping$Nb_ttx
data_phenotyping$p_larvae <- data_phenotyping$Larvae / data_phenotyping$Nb_ttx
data_phenotyping_all <- data_phenotyping
data_phenotyping <- data_phenotyping[data_phenotyping$Week=="5-week",]1.4 Merge dataset He
#Add heterozygosity
data <- merge(x = data, y = data_he, by = "Population", all.x=TRUE)
#Add heterozygosity
#To add population extinct
data_phenotyping_withextinct <- merge(x = data_phenotyping, y = data_he, by = "Population", all.y=TRUE)
data_phenotyping <- merge(x = data_phenotyping, y = data_he, by = "Population", all.x=TRUE)1.5 Subset dataset generation
#Data beginning
data_G2 <- data[data$Generation_Eggs==2,]
#Data end
data_G6 <- data[data$Generation_Eggs==6,]1.6 Merge dataset phenotyping
## Merge dataset phenotyping
#merge dataset
temp_Start <- data_G2[,c("Block","Population",
"Genetic_Diversity","Lambda", "He_remain", "He_remain_end" )]
temp_Start$ID_Rep <- 1
temp_Start$Phenotyping <- "Initial"
temp_end <- data_phenotyping[,c("Block","Population",
"Genetic_Diversity","Lambda","ID_Rep", "He_remain", "He_remain_end" )]
temp_end$Phenotyping <- "Final"
data_evolution_Lambda <- rbind(temp_Start,temp_end)
data_evolution_Lambda$Phenotyping <- factor(data_evolution_Lambda$Phenotyping,
levels = c("Initial", "Final"))2 PLOT
2.1 Population size
PLOT_Pop_size<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Population,
colour = Genetic_Diversity)) +
geom_point(position = position_dodge(0.4), size = 0.8, alpha = 0.7) +
geom_line(position = position_dodge(0.4), size = 0.2, alpha = 0.6) +
ylab("Population size") +
labs(color="Genetic diversity") +
xlab("Generation") +
theme_LO
PLOT_Pop_size
# cowplot::save_plot(file = here::here("figures","OldPlot_PopulationSize.pdf"),PLOT_Pop_size, base_height = 10/cm(1),
# base_width = 15/cm(1), dpi = 610)
SUM_Popsize <- Rmisc::summarySE(data,
measurevar = c("Nb_adults"),
groupvars = c("Generation_Eggs",
"Genetic_Diversity"),
na.rm=TRUE)
SUM_Popsize## Generation_Eggs Genetic_Diversity N Nb_adults sd se ci
## 1 1 High 27 100.00000 0.00000 0.000000 0.00000
## 2 1 Medium 23 100.00000 0.00000 0.000000 0.00000
## 3 1 Low 34 100.00000 0.00000 0.000000 0.00000
## 4 2 High 27 344.29630 73.77294 14.197610 29.18360
## 5 2 Medium 23 282.04348 100.75735 21.009360 43.57075
## 6 2 Low 34 282.61765 78.53006 13.467794 27.40043
## 7 3 High 27 151.85185 80.96899 15.582489 32.03027
## 8 3 Medium 23 118.73913 88.54491 18.462891 38.28969
## 9 3 Low 34 104.61765 85.13341 14.600260 29.70445
## 10 4 High 26 114.00000 80.20224 15.728954 32.39439
## 11 4 Medium 23 62.65217 64.50483 13.450188 27.89398
## 12 4 Low 34 46.38235 39.63241 6.796903 13.82840
## 13 5 High 27 103.62963 51.56486 9.923661 20.39838
## 14 5 Medium 20 61.50000 50.59904 11.314290 23.68108
## 15 5 Low 31 51.51613 46.13377 8.285870 16.92200
## 16 6 High 27 147.66667 42.60282 8.198916 16.85311
## 17 6 Medium 19 69.73684 46.02396 10.558620 22.18284
## 18 6 Low 28 77.03571 46.14428 8.720450 17.89289# SUM_Popsize <- Rmisc::summarySE(data[data$Extinction_finalstatus=="no",],
# measurevar = c("Nb_adults"),
# groupvars = c("Generation_Eggs",
# "Genetic_Diversity"),
# na.rm=TRUE)
# SUM_Popsize
PLOT_Pop_size_average <- PLOT_Pop_size +
geom_point(data = SUM_Popsize, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Genetic_Diversity,
colour = Genetic_Diversity),
size = 5, position = position_dodge(0.4)) +
geom_errorbar(data = SUM_Popsize, aes(x = factor(Generation_Eggs),
ymin = Nb_adults-ci, ymax = Nb_adults+ci,
group = Genetic_Diversity,
colour = Genetic_Diversity),
size = 1, width=.2, position = position_dodge(0.4)) +
geom_line(data = SUM_Popsize, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Genetic_Diversity,
colour = Genetic_Diversity),
linetype = "dashed", size = 1.5, position = position_dodge(0.4))
PLOT_Pop_size_average
# cowplot::save_plot(file = here::here("figures","PopulationSize_average.pdf"), PLOT_Pop_size_average, base_height = 10/cm(1), base_width = 15/cm(1), dpi = 610)2.2 Extinction
## Extinction yes no
PLOT_Extinction<-ggplot2::ggplot(data, aes(x = factor(Generation_Eggs),
y = Nb_adults,
group = Population,
colour = Extinction_finalstatus)) +
facet_wrap(~Genetic_Diversity, ncol=3) +
#geom_point(position = position_dodge(0.4), size = 0.4, alpha = 0.7) +
geom_line(position = position_dodge(0.1), size = 0.4, alpha = 0.8) +
ylab("Population size") +
scale_color_manual(values = c("black", "red")) +
xlab("Generation") +
theme(legend.position = "none") +
theme_LO
PLOT_Extinction
# #
# cowplot::save_plot(here::here("figures","Extinction.pdf"), PLOT_Extinction, base_height = 6/cm(1),
# base_width = 20/cm(1), dpi = 610)
#
# #2.3 Growth rate G2
tapply(data_G2$Lambda,data_G2$Genetic_Diversity,mean)## High Medium Low
## 3.442963 2.820435 2.826176PLOT_Growth<-ggplot2::ggplot(data_G2, aes(x = Genetic_Diversity, y = Lambda,
colour = Genetic_Diversity)) +
geom_boxplot(aes(group = Genetic_Diversity),outlier.colour = NA) +
geom_jitter(width = 0.25, size = 1, alpha = 0.8) +
ylab("Growth rate") +
theme(legend.position = "none") +
xlab("Genetic diversity") +
theme_LO
PLOT_Growth
#
# cowplot::save_plot(here::here("figures","OldPlot_Growth_Start.pdf"), PLOT_Growth, base_height = 10/cm(1),
# base_width = 8/cm(1), dpi = 610)2.4 Growth rate evolution
SUM_MEAN_start_end <- Rmisc::summarySE(data_evolution_Lambda, measurevar = c("Lambda"),
groupvars = c("Genetic_Diversity","Phenotyping"),
na.rm=TRUE)
SUM_POP_start_end <- Rmisc::summarySE(data_evolution_Lambda, measurevar = c("Lambda"),
groupvars = c("Genetic_Diversity","Phenotyping",
"Population"),
na.rm=TRUE)
PLOT_Growth_startend<-ggplot(SUM_POP_start_end, aes(x = Phenotyping,
y = Lambda,
group = Population,
colour = Genetic_Diversity)) +
geom_point(position = position_dodge(0.4), size = 0.8, alpha = 0.7) +
geom_line(position = position_dodge(0.4), size = 0.2, alpha = 0.6) +
ylab("Growth rate") +
xlab("Phenotyping") +
theme_LO
PLOT_Growth_startend
PLOT_Growth_startend_average <- PLOT_Growth_startend +
geom_point(data = SUM_MEAN_start_end, aes(x = Phenotyping,
y = Lambda,
group = Genetic_Diversity,
colour = Genetic_Diversity),
size = 5, position = position_dodge(0.4)) +
geom_errorbar(data = SUM_MEAN_start_end, aes(x = Phenotyping,
ymin = Lambda-ci, ymax = Lambda+ci,
group = Genetic_Diversity,
colour = Genetic_Diversity),
size = 1, width=.2, position = position_dodge(0.4)) +
geom_line(data = SUM_MEAN_start_end, aes(x = Phenotyping,
y = Lambda,
group = Genetic_Diversity,
colour = Genetic_Diversity),
linetype = "dashed", size = 1.5, position = position_dodge(0.4))
PLOT_Growth_startend_average
#
# cowplot::save_plot(here::here("figures","OldPlot_Lambda_Evolution.pdf"),
# PLOT_Growth_startend_average, base_height = 10/cm(1),
# base_width = 15/cm(1), dpi = 610)2.5 Life stage proportion
SUM_Prop_adults<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_adults"),
groupvars = c("Genetic_Diversity","Week"),
na.rm = TRUE)
SUM_Prop_adults## Genetic_Diversity Week N p_adults sd se ci
## 1 High 4-week 50 0.6508098 0.17776618 0.025139934 0.05052059
## 2 High 5-week 105 0.9632040 0.05259208 0.005132461 0.01017786
## 3 Medium 4-week 24 0.7121270 0.15722597 0.032093616 0.06639070
## 4 Medium 5-week 49 0.9563136 0.06160975 0.008801393 0.01769639
## 5 Low 4-week 30 0.6568606 0.19029917 0.034743716 0.07105888
## 6 Low 5-week 59 0.9582021 0.07845154 0.010213520 0.02044458SUM_Prop_pupae<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_pupae"),
groupvars = c("Genetic_Diversity","Week"),
na.rm = TRUE)
SUM_Prop_pupae## Genetic_Diversity Week N p_pupae sd se ci
## 1 High 4-week 50 0.28825555 0.15629878 0.022103985 0.044419621
## 2 High 5-week 105 0.01458769 0.02407445 0.002349426 0.004658999
## 3 Medium 4-week 24 0.23019148 0.14506173 0.029610601 0.061254195
## 4 Medium 5-week 49 0.02283167 0.03485130 0.004978757 0.010010462
## 5 Low 4-week 30 0.29380227 0.17562943 0.032065401 0.065581108
## 6 Low 5-week 59 0.01675308 0.02459640 0.003202178 0.006409856SUM_Prop_larvae<- Rmisc::summarySE(data_phenotyping_all,
measurevar = c("p_larvae"),
groupvars = c("Genetic_Diversity","Week"),
na.rm = TRUE)
SUM_Prop_larvae## Genetic_Diversity Week N p_larvae sd se ci
## 1 High 4-week 50 0.06093466 0.04274850 0.006045551 0.012148989
## 2 High 5-week 105 0.02220828 0.04324043 0.004219833 0.008368088
## 3 Medium 4-week 24 0.05768150 0.05999020 0.012245449 0.025331640
## 4 Medium 5-week 49 0.02085474 0.04462923 0.006375605 0.012819012
## 5 Low 4-week 30 0.04933713 0.08198954 0.014969174 0.030615398
## 6 Low 5-week 59 0.02504479 0.06691319 0.008711355 0.017437671SUM_Prop <- data.frame(SUM_Prop_adults[,1:4],
p_pupae=SUM_Prop_pupae$p_pupae,
p_larvae=SUM_Prop_larvae$p_larvae)
SUM_Prop<-tidyr::gather(SUM_Prop,"Stage", "Proportion",4:6,-"Genetic_Diversity")
SUM_Prop## Genetic_Diversity Week N Stage Proportion
## 1 High 4-week 50 p_adults 0.65080979
## 2 High 5-week 105 p_adults 0.96320403
## 3 Medium 4-week 24 p_adults 0.71212702
## 4 Medium 5-week 49 p_adults 0.95631359
## 5 Low 4-week 30 p_adults 0.65686059
## 6 Low 5-week 59 p_adults 0.95820213
## 7 High 4-week 50 p_pupae 0.28825555
## 8 High 5-week 105 p_pupae 0.01458769
## 9 Medium 4-week 24 p_pupae 0.23019148
## 10 Medium 5-week 49 p_pupae 0.02283167
## 11 Low 4-week 30 p_pupae 0.29380227
## 12 Low 5-week 59 p_pupae 0.01675308
## 13 High 4-week 50 p_larvae 0.06093466
## 14 High 5-week 105 p_larvae 0.02220828
## 15 Medium 4-week 24 p_larvae 0.05768150
## 16 Medium 5-week 49 p_larvae 0.02085474
## 17 Low 4-week 30 p_larvae 0.04933713
## 18 Low 5-week 59 p_larvae 0.02504479SUM_Prop$Stage <- factor(SUM_Prop$Stage, levels = c("p_adults", "p_pupae", "p_larvae"))
PLOT_prop <- ggplot2::ggplot(data = SUM_Prop, aes(x = Genetic_Diversity,
y = Proportion, fill = Stage)) +
facet_wrap(~ Week, ncol=2) +
geom_bar(stat="identity") +
xlab("Level of genetic diversity") +
labs(color="Stage of individuals") +
ylab("Proportion of each stage") +
scale_fill_manual(values= c("#4EABB9", "#E95C2B", "#CBA938"),
breaks = c("p_adults", "p_pupae", "p_larvae"),
labels = c( "Adult", "Pupae","Larvae")) +
theme_LO
PLOT_prop
#
# cowplot::save_plot(here::here("figures","Life_stage_proportion.pdf"), PLOT_prop,
# base_height = 8/cm(1), base_width = 16/cm(1), dpi = 610)2.6 Remaining heterozygosity
# All but with: He beginning for both and with pseudoreplication for final
ggplot2::ggplot(data_evolution_Lambda, aes(x = He_remain, y = Lambda,
colour = Genetic_Diversity)) +
facet_wrap(~ Phenotyping, ncol = 1, scales = "free")+
geom_point(size = 2, alpha = 0.8) +
ylab("Growth rate") +
theme(legend.position = "none") +
xlab("Remaining heterozygosity") +
theme_LO
PLOT_He_Initial <- ggplot2::ggplot(data_G2, aes(x = He_remain, y = Lambda,
colour = Genetic_Diversity)) +
geom_point(size = 2, alpha = 0.8) +
ylab("Growth rate") +
theme(legend.position = "none") +
xlab("Remaining heterozygosity") +
theme_LO
PLOT_He_Initial
PLOT_He_End <- ggplot2::ggplot(data_evolution_Lambda[data_evolution_Lambda$Phenotyping=="Final",],
aes(x = He_remain_end, y = Lambda,
colour = Genetic_Diversity)) +
geom_point(size = 2, alpha = 0.8) +
ylab("Growth rate") +
theme(legend.position = "none") +
xlab("Remaining heterozygosity") +
theme_LO
PLOT_He_End
# Remove pseudoreplication: average of each population
SUM_POP_He_Mean<- Rmisc::summarySE(data_evolution_Lambda[data_evolution_Lambda$Phenotyping=="Final",],
measurevar = c("Lambda"),
groupvars = c("Genetic_Diversity","He_remain_end",
"Population"),
na.rm=TRUE)
PLOT_He_Final <- ggplot2::ggplot(SUM_POP_He_Mean, aes(x = He_remain_end, y = Lambda,
colour = Genetic_Diversity)) +
geom_point(size = 2, alpha = 0.8) +
ylab("Growth rate") +
theme(legend.position = "none") +
xlab("Remaining heterozygosity") +
theme_LO
PLOT_He_Final
# ALL WITHOUT PSEUDO
SUM_POP_He_ALL<- Rmisc::summarySE(data_evolution_Lambda,
measurevar = c("Lambda"),
groupvars = c("Genetic_Diversity","Phenotyping", "He_remain", "He_remain_end",
"Population"),
na.rm=TRUE)
SUM_POP_He_ALL$HE <- ifelse(SUM_POP_He_ALL$Phenotyping=="Initial",SUM_POP_He_ALL$He_remain,SUM_POP_He_ALL$He_remain_end)
SUM_POP_He_ALL <- SUM_POP_He_ALL[!is.na(SUM_POP_He_ALL$He_remain_end),]
SUM_POP_He_ALL <- SUM_POP_He_ALL[SUM_POP_He_ALL$HE!="-Inf",]
facet_names <- c('Initial'="Before rescue experiment",'Final'="After rescue experiment")
PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda,
colour = Genetic_Diversity,
shape = Phenotyping)) +
#facet_wrap(~ Phenotyping, ncol=1) +
geom_point(size = 3, alpha = 0.8) +
scale_shape_manual(values = c(1, 19)) +
ylab("Growth rate") +
xlab("Remaining heterozygosity") +
theme_LO
PLOT_He_ALL
PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda,
colour = Genetic_Diversity,
shape = Phenotyping)) +
facet_wrap(~ Phenotyping, ncol=1) +
geom_point(size = 3, alpha = 0.8) +
scale_shape_manual(values = c(1, 19)) +
ylab("Growth rate") +
xlab("Remaining heterozygosity") +
theme_LO
PLOT_He_ALL
PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL, aes(x = HE, y = Lambda,
colour = Genetic_Diversity)) +
facet_wrap(~ Phenotyping, ncol=1) +
geom_point(size = 3, alpha = 0.8) +
ylab("Growth rate") +
xlab("Remaining heterozygosity") +
theme_LO
PLOT_He_ALL
# cowplot::save_plot(here::here("figures","Fitness_He_initial_final.pdf"), PLOT_He_ALL,
# base_height = 20/cm(1), base_width = 16/cm(1), dpi = 610)
# PLOT_He_ALL <- ggplot2::ggplot(SUM_POP_He_ALL[SUM_POP_He_ALL$Phenotyping=="Final",], aes(x = HE, y = Lambda,
# colour = Genetic_Diversity,
# shape = Phenotyping)) +
# geom_point(size = 3, alpha = 0.8) +
# scale_shape_manual(values = c(19)) +
# ylab("Growth rate") +
# xlab("Remaining heterozygosity") +
# theme_LO
# PLOT_He_ALL
########################
###Plot final phenotyping with extinct
########################
data_phenotyping_withextinct## Population Week Block ID_Rep Divsersity Population.1 Box Initial.N
## 1 High1 5-week Block3 52 High 1 1 30
## 2 High1 5-week Block3 53 High 1 2 30
## 3 High13 5-week Block4 130 High 13 2 30
## 4 High13 5-week Block4 131 High 13 3 30
## 5 High13 5-week Block4 129 High 13 1 30
## 6 High14 5-week Block4 132 High 14 1 30
## 7 High14 5-week Block4 133 High 14 2 30
## 8 High14 5-week Block4 134 High 14 3 30
## 9 High14 5-week Block4 135 High 14 4 30
## 10 High15 5-week Block4 136 High 15 1 30
## 11 High15 5-week Block4 137 High 15 2 30
## 12 High16 5-week Block4 138 High 16 1 30
## 13 High16 5-week Block4 139 High 16 2 30
## 14 High16 5-week Block4 140 High 16 3 30
## 15 High16 5-week Block4 141 High 16 4 30
## 16 High16 5-week Block4 142 High 16 5 30
## 17 High17 5-week Block4 143 High 17 1 30
## 18 High18 5-week Block4 145 High 18 2 30
## 19 High18 5-week Block4 146 High 18 3 30
## 20 High18 5-week Block4 144 High 18 1 30
## 21 High19 5-week Block4 147 High 19 1 30
## 22 High19 5-week Block4 149 High 19 3 30
## 23 High19 5-week Block4 148 High 19 2 30
## 24 High2 5-week Block3 54 High 2 1 30
## 25 High2 5-week Block3 55 High 2 2 30
## 26 High20 5-week Block4 150 High 20 1 30
## 27 High20 5-week Block4 151 High 20 2 30
## 28 High20 5-week Block4 152 High 20 3 30
## 29 High20 5-week Block4 153 High 20 4 30
## 30 High20 5-week Block4 154 High 20 5 30
## 31 High21 5-week Block4 155 High 21 1 30
## 32 High21 5-week Block4 157 High 21 3 30
## 33 High21 5-week Block4 156 High 21 2 30
## 34 High22 5-week Block4 158 High 22 1 30
## 35 High22 5-week Block4 159 High 22 2 30
## 36 High23 5-week Block4 160 High 23 1 30
## 37 High23 5-week Block4 162 High 23 3 30
## 38 High23 5-week Block4 163 High 23 4 30
## 39 High23 5-week Block4 161 High 23 2 30
## 40 High24 5-week Block4 164 High 24 1 30
## 41 High24 5-week Block4 165 High 24 2 30
## 42 High25 5-week Block5 185 High 25 3 30
## 43 High25 5-week Block5 184 High 25 2 30
## 44 High25 5-week Block5 183 High 25 1 30
## 45 High25 5-week Block5 186 High 25 4 30
## 46 High25 5-week Block5 187 High 25 5 30
## 47 High26 5-week Block5 188 High 26 1 30
## 48 High27 5-week Block5 190 High 27 2 30
## 49 High27 5-week Block5 189 High 27 1 30
## 50 High27 5-week Block5 191 High 27 3 30
## 51 High27 5-week Block5 192 High 27 4 30
## 52 High28 5-week Block5 198 High 28 6 30
## 53 High28 5-week Block5 193 High 28 1 30
## 54 High28 5-week Block5 196 High 28 4 30
## 55 High28 5-week Block5 195 High 28 3 30
## 56 High28 5-week Block5 194 High 28 2 30
## 57 High28 5-week Block5 197 High 28 5 30
## 58 High3 5-week Block3 56 High 3 1 30
## 59 High3 5-week Block3 57 High 3 2 30
## 60 High3 5-week Block3 58 High 3 3 30
## 61 High32 5-week Block5 203 High 32 1 30
## 62 High32 5-week Block5 204 High 32 2 30
## 63 High33 5-week Block5 205 High 33 1 30
## 64 High33 5-week Block5 206 High 33 2 30
## 65 High34 5-week Block5 210 High 34 4 30
## 66 High34 5-week Block5 211 High 34 5 30
## 67 High34 5-week Block5 209 High 34 3 30
## 68 High34 5-week Block5 207 High 34 1 30
## 69 High34 5-week Block5 208 High 34 2 30
## 70 High35 5-week Block5 215 High 35 4 30
## 71 High35 5-week Block5 212 High 35 1 30
## 72 High35 5-week Block5 213 High 35 2 30
## 73 High35 5-week Block5 214 High 35 3 30
## 74 High35 5-week Block5 217 High 35 6 30
## 75 High35 5-week Block5 216 High 35 5 30
## 76 High4 5-week Block3 60 High 4 2 30
## 77 High4 5-week Block3 61 High 4 3 30
## 78 High4 5-week Block3 62 High 4 4 30
## 79 High4 5-week Block3 59 High 4 1 30
## 80 High4 5-week Block3 63 High 4 5 30
## 81 High5 5-week Block3 64 High 5 1 30
## 82 High5 5-week Block3 65 High 5 2 30
## 83 High5 5-week Block3 66 High 5 3 30
## 84 High5 5-week Block3 67 High 5 4 30
## 85 High5 5-week Block3 68 High 5 5 30
## 86 High5 5-week Block3 69 High 5 6 30
## 87 High5 5-week Block3 70 High 5 7 30
## 88 High6 5-week Block3 72 High 6 2 30
## 89 High6 5-week Block3 73 High 6 3 30
## 90 High6 5-week Block3 71 High 6 1 30
## 91 High7 5-week Block3 74 High 7 1 30
## 92 High7 5-week Block3 75 High 7 2 30
## 93 High7 5-week Block3 77 High 7 4 30
## 94 High7 5-week Block3 78 High 7 5 30
## 95 High7 5-week Block3 79 High 7 6 30
## 96 High7 5-week Block3 76 High 7 3 30
## 97 Low1 5-week Block3 5 Low 1 5 30
## 98 Low1 5-week Block3 1 Low 1 1 30
## 99 Low1 5-week Block3 2 Low 1 2 30
## 100 Low1 5-week Block3 3 Low 1 3 30
## 101 Low1 5-week Block3 4 Low 1 4 30
## 102 Low10 5-week Block3 29 Low 10 1 30
## 103 Low11 5-week Block3 30 Low 11 1 30
## 104 Low12 5-week Block3 32 Low 12 2 30
## 105 Low12 5-week Block3 31 Low 12 1 30
## 106 Low13 5-week Block4 85 Low 13 1 30
## 107 Low14 5-week Block4 87 Low 14 2 30
## 108 Low14 5-week Block4 88 Low 14 3 30
## 109 Low14 5-week Block4 86 Low 14 1 30
## 110 Low15 5-week Block4 90 Low 15 2 30
## 111 Low15 5-week Block4 89 Low 15 1 30
## 112 Low16 5-week Block4 91 Low 16 1 30
## 113 Low17 5-week Block4 93 Low 17 2 30
## 114 Low17 5-week Block4 92 Low 17 1 30
## 115 Low18 5-week Block4 97 Low 18 4 30
## 116 Low18 5-week Block4 94 Low 18 1 30
## 117 Low18 5-week Block4 95 Low 18 2 30
## 118 Low18 5-week Block4 96 Low 18 3 30
## 119 Low19 5-week Block4 98 Low 19 1 30
## 120 Low2 5-week Block3 7 Low 2 2 30
## 121 Low2 5-week Block3 8 Low 2 3 30
## 122 Low2 5-week Block3 9 Low 2 4 30
## 123 Low2 5-week Block3 6 Low 2 1 30
## 124 Low20 5-week Block4 99 Low 20 1 30
## 125 Low21 5-week Block4 100 Low 21 1 30
## 126 Low22 5-week Block4 101 Low 22 1 30
## 127 Low24 5-week Block4 103 Low 24 1 30
## 128 Low25 5-week Block4 105 Low 25 2 30
## 129 Low25 5-week Block4 104 Low 25 1 30
## 130 Low26 5-week Block5 166 Low 26 1 30
## 131 Low29 5-week Block5 167 Low 29 1 30
## 132 Low3 5-week Block3 11 Low 3 2 30
## 133 Low3 5-week Block3 12 Low 3 3 30
## 134 Low3 5-week Block3 10 Low 3 1 30
## 135 Low33 5-week Block5 168 Low 33 1 30
## 136 Low34 5-week Block5 171 Low 34 3 30
## 137 Low34 5-week Block5 172 Low 34 4 30
## 138 Low34 5-week Block5 170 Low 34 2 30
## 139 Low34 5-week Block5 169 Low 34 1 30
## 140 Low35 5-week Block5 173 Low 35 1 30
## 141 Low4 5-week Block3 13 Low 4 1 30
## 142 Low4 5-week Block3 14 Low 4 2 30
## 143 Low4 5-week Block3 15 Low 4 3 30
## 144 Low4 5-week Block3 16 Low 4 4 30
## 145 Low4 5-week Block3 17 Low 4 5 30
## 146 Low5 5-week Block3 19 Low 5 2 30
## 147 Low5 5-week Block3 20 Low 5 3 30
## 148 Low5 5-week Block3 18 Low 5 1 30
## 149 Low8 5-week Block3 21 Low 8 1 30
## 150 Low8 5-week Block3 22 Low 8 2 30
## 151 Low8 5-week Block3 24 Low 8 4 30
## 152 Low8 5-week Block3 23 Low 8 3 30
## 153 Low9 5-week Block3 25 Low 9 1 30
## 154 Low9 5-week Block3 26 Low 9 2 30
## 155 Low9 5-week Block3 27 Low 9 3 30
## 156 Low9 5-week Block3 28 Low 9 4 30
## 157 Med1 5-week Block3 33 Med 1 1 30
## 158 Med10 5-week Block4 110 Med 10 1 30
## 159 Med11 5-week Block4 112 Med 11 2 30
## 160 Med11 5-week Block4 113 Med 11 3 30
## 161 Med11 5-week Block4 111 Med 11 1 30
## 162 Med11 5-week Block4 114 Med 11 4 30
## 163 Med12 5-week Block4 115 Med 12 1 30
## 164 Med13 5-week Block4 118 Med 13 3 30
## 165 Med13 5-week Block4 116 Med 13 1 30
## 166 Med13 5-week Block4 117 Med 13 2 30
## 167 Med15 5-week Block4 119 Med 15 1 30
## 168 Med15 5-week Block4 122 Med 15 4 30
## 169 Med15 5-week Block4 123 Med 15 5 30
## 170 Med15 5-week Block4 120 Med 15 2 30
## 171 Med15 5-week Block4 121 Med 15 3 30
## 172 Med18 5-week Block5 174 Med 18 1 30
## 173 Med18 5-week Block5 176 Med 18 3 30
## 174 Med18 5-week Block5 177 Med 18 4 30
## 175 Med18 5-week Block5 175 Med 18 2 30
## 176 Med19 5-week Block5 178 Med 19 1 30
## 177 Med2 5-week Block3 34 Med 2 1 30
## 178 Med2 5-week Block3 35 Med 2 2 30
## 179 Med2 5-week Block3 36 Med 2 3 30
## 180 Med21 5-week Block5 179 Med 21 1 30
## 181 Med23 5-week Block5 180 Med 23 1 30
## 182 Med24 5-week Block5 181 Med 24 1 30
## 183 Med24 5-week Block5 182 Med 24 2 30
## 184 Med3 5-week Block3 37 Med 3 1 30
## 185 Med3 5-week Block3 38 Med 3 2 30
## 186 Med3 5-week Block3 39 Med 3 3 30
## 187 Med4 5-week Block3 40 Med 4 1 30
## 188 Med5 5-week Block3 41 Med 5 1 30
## 189 Med6 5-week Block3 44 Med 6 3 30
## 190 Med6 5-week Block3 42 Med 6 1 30
## 191 Med6 5-week Block3 43 Med 6 2 30
## 192 Med7 5-week Block3 46 Med 7 2 30
## 193 Med7 5-week Block3 45 Med 7 1 30
## 194 Med8 5-week Block3 47 Med 8 1 30
## 195 Med8 5-week Block3 48 Med 8 2 30
## 196 Med8 5-week Block3 49 Med 8 3 30
## 197 Med8 5-week Block3 51 Med 8 5 30
## 198 Med8 5-week Block3 50 Med 8 4 30
## 199 Med9 5-week Block4 106 Med 9 1 30
## 200 Med9 5-week Block4 109 Med 9 4 30
## 201 Med9 5-week Block4 107 Med 9 2 30
## 202 Med9 5-week Block4 108 Med 9 3 30
## 203 Low27 <NA> <NA> <NA> <NA> NA NA NA
## 204 Low28 <NA> <NA> <NA> <NA> NA NA NA
## 205 Low30 <NA> <NA> <NA> <NA> NA NA NA
## 206 Low31 <NA> <NA> <NA> <NA> NA NA NA
## 207 Low32 <NA> <NA> <NA> <NA> NA NA NA
## 208 Low36 <NA> <NA> <NA> <NA> NA NA NA
## 209 Low7 <NA> <NA> <NA> <NA> NA NA NA
## 210 Med14 <NA> <NA> <NA> <NA> NA NA NA
## 211 Med17 <NA> <NA> <NA> <NA> NA NA NA
## 212 Med20 <NA> <NA> <NA> <NA> NA NA NA
## 213 Med22 <NA> <NA> <NA> <NA> NA NA NA
## Start End Larvae Pupae Adults Lambda Genetic_Diversity Nb_ttx
## 1 7/8/21 8/12/21 3 9 97 3.23333333 High 109
## 2 7/8/21 8/12/21 1 1 8 0.26666667 High 10
## 3 7/15/21 8/19/21 1 2 104 3.46666667 High 107
## 4 7/15/21 8/19/21 2 1 84 2.80000000 High 87
## 5 7/15/21 8/19/21 2 0 84 2.80000000 High 86
## 6 7/15/21 8/19/21 1 0 83 2.76666667 High 84
## 7 7/15/21 8/19/21 2 0 48 1.60000000 High 50
## 8 7/15/21 8/19/21 0 0 62 2.06666667 High 62
## 9 7/15/21 8/19/21 1 0 43 1.43333333 High 44
## 10 7/15/21 8/19/21 4 0 84 2.80000000 High 88
## 11 7/15/21 8/19/21 8 0 57 1.90000000 High 65
## 12 7/15/21 8/19/21 1 1 110 3.66666667 High 112
## 13 7/15/21 8/19/21 1 0 45 1.50000000 High 46
## 14 7/15/21 8/19/21 0 0 59 1.96666667 High 59
## 15 7/15/21 8/19/21 3 0 54 1.80000000 High 57
## 16 7/15/21 8/19/21 1 2 56 1.86666667 High 59
## 17 7/15/21 8/19/21 0 3 165 5.50000000 High 168
## 18 7/15/21 8/19/21 3 2 50 1.66666667 High 55
## 19 7/15/21 8/19/21 0 0 88 2.93333333 High 88
## 20 7/15/21 8/19/21 0 5 75 2.50000000 High 80
## 21 7/15/21 8/19/21 3 6 66 2.20000000 High 75
## 22 7/15/21 8/19/21 0 1 93 3.10000000 High 94
## 23 7/15/21 8/19/21 1 1 112 3.73333333 High 114
## 24 7/8/21 8/12/21 4 0 165 5.50000000 High 169
## 25 7/8/21 8/12/21 0 1 73 2.43333333 High 74
## 26 7/15/21 8/19/21 5 0 77 2.56666667 High 82
## 27 7/15/21 8/19/21 1 1 84 2.80000000 High 86
## 28 7/15/21 8/19/21 0 0 101 3.36666667 High 101
## 29 7/15/21 8/19/21 2 1 48 1.60000000 High 51
## 30 7/15/21 8/19/21 0 0 80 2.66666667 High 80
## 31 7/15/21 8/19/21 0 0 92 3.06666667 High 92
## 32 7/15/21 8/19/21 1 4 86 2.86666667 High 91
## 33 7/15/21 8/19/21 2 2 121 4.03333333 High 125
## 34 7/15/21 8/19/21 2 4 135 4.50000000 High 141
## 35 7/15/21 8/19/21 0 0 83 2.76666667 High 83
## 36 7/15/21 8/19/21 1 0 61 2.03333333 High 62
## 37 7/15/21 8/19/21 0 0 86 2.86666667 High 86
## 38 7/15/21 8/19/21 0 0 47 1.56666667 High 47
## 39 7/15/21 8/19/21 0 0 66 2.20000000 High 66
## 40 7/15/21 8/19/21 2 5 88 2.93333333 High 95
## 41 7/15/21 8/19/21 12 3 105 3.50000000 High 120
## 42 7/22/21 8/26/21 2 1 98 3.26666667 High 101
## 43 7/22/21 8/26/21 1 5 71 2.36666667 High 77
## 44 7/22/21 8/26/21 2 1 59 1.96666667 High 62
## 45 7/22/21 8/26/21 4 2 86 2.86666667 High 92
## 46 7/22/21 8/26/21 1 0 91 3.03333333 High 92
## 47 7/22/21 8/26/21 0 0 79 2.63333333 High 79
## 48 7/22/21 8/26/21 4 3 71 2.36666667 High 78
## 49 7/22/21 8/26/21 37 4 73 2.43333333 High 114
## 50 7/22/21 8/26/21 1 2 73 2.43333333 High 76
## 51 7/22/21 8/26/21 2 7 83 2.76666667 High 92
## 52 7/22/21 8/26/21 0 0 23 0.76666667 High 23
## 53 7/22/21 8/26/21 1 1 103 3.43333333 High 105
## 54 7/22/21 8/26/21 1 3 65 2.16666667 High 69
## 55 7/22/21 8/26/21 1 0 75 2.50000000 High 76
## 56 7/22/21 8/26/21 1 0 129 4.30000000 High 130
## 57 7/22/21 8/26/21 0 0 56 1.86666667 High 56
## 58 7/8/21 8/12/21 1 0 124 4.13333333 High 125
## 59 7/8/21 8/12/21 16 0 72 2.40000000 High 88
## 60 7/8/21 8/12/21 0 0 35 1.16666667 High 35
## 61 7/22/21 8/26/21 3 2 136 4.53333333 High 141
## 62 7/22/21 8/26/21 2 0 143 4.76666667 High 145
## 63 7/22/21 8/26/21 1 2 61 2.03333333 High 64
## 64 7/22/21 8/26/21 2 0 36 1.20000000 High 38
## 65 7/22/21 8/26/21 0 0 116 3.86666667 High 116
## 66 7/22/21 8/26/21 0 0 36 1.20000000 High 36
## 67 7/22/21 8/26/21 0 1 70 2.33333333 High 71
## 68 7/22/21 8/26/21 1 1 101 3.36666667 High 103
## 69 7/22/21 8/26/21 0 0 82 2.73333333 High 82
## 70 7/22/21 8/26/21 1 0 101 3.36666667 High 102
## 71 7/22/21 8/26/21 0 0 107 3.56666667 High 107
## 72 7/22/21 8/26/21 0 0 49 1.63333333 High 49
## 73 7/22/21 8/26/21 1 0 75 2.50000000 High 76
## 74 7/22/21 8/26/21 0 0 36 1.20000000 High 36
## 75 7/22/21 8/26/21 0 0 48 1.60000000 High 48
## 76 7/8/21 8/12/21 0 0 186 6.20000000 High 186
## 77 7/8/21 8/12/21 0 0 107 3.56666667 High 107
## 78 7/8/21 8/12/21 0 0 73 2.43333333 High 73
## 79 7/8/21 8/12/21 0 0 145 4.83333333 High 145
## 80 7/8/21 8/12/21 0 0 35 1.16666667 High 35
## 81 7/8/21 8/12/21 0 1 44 1.46666667 High 45
## 82 7/8/21 8/12/21 0 0 93 3.10000000 High 93
## 83 7/8/21 8/12/21 0 0 84 2.80000000 High 84
## 84 7/8/21 8/12/21 1 0 81 2.70000000 High 82
## 85 7/8/21 8/12/21 0 0 66 2.20000000 High 66
## 86 7/8/21 8/12/21 0 0 51 1.70000000 High 51
## 87 7/8/21 8/12/21 0 1 36 1.20000000 High 37
## 88 7/8/21 8/12/21 1 0 61 2.03333333 High 62
## 89 7/8/21 8/12/21 1 0 17 0.56666667 High 18
## 90 7/8/21 8/12/21 2 1 66 2.20000000 High 69
## 91 7/8/21 8/12/21 0 1 101 3.36666667 High 102
## 92 7/8/21 8/12/21 3 2 97 3.23333333 High 102
## 93 7/8/21 8/12/21 0 8 87 2.90000000 High 95
## 94 7/8/21 8/12/21 0 0 113 3.76666667 High 113
## 95 7/8/21 8/12/21 0 4 46 1.53333333 High 50
## 96 7/8/21 8/12/21 0 0 82 2.73333333 High 82
## 97 7/8/21 8/12/21 0 0 23 0.76666667 Low 23
## 98 7/8/21 8/12/21 1 0 44 1.46666667 Low 45
## 99 7/8/21 8/12/21 0 0 42 1.40000000 Low 42
## 100 7/8/21 8/12/21 0 0 49 1.63333333 Low 49
## 101 7/8/21 8/12/21 1 1 57 1.90000000 Low 59
## 102 7/8/21 8/12/21 0 2 56 1.86666667 Low 58
## 103 7/8/21 8/12/21 0 0 17 0.56666667 Low 17
## 104 7/8/21 8/12/21 1 0 57 1.90000000 Low 58
## 105 7/8/21 8/12/21 2 1 117 3.90000000 Low 120
## 106 7/15/21 8/19/21 6 1 77 2.56666667 Low 84
## 107 7/15/21 8/19/21 0 1 96 3.20000000 Low 97
## 108 7/15/21 8/19/21 2 2 91 3.03333333 Low 95
## 109 7/15/21 8/19/21 2 1 102 3.40000000 Low 105
## 110 7/15/21 8/19/21 1 3 72 2.40000000 Low 76
## 111 7/15/21 8/19/21 3 1 51 1.70000000 Low 55
## 112 7/15/21 8/19/21 0 0 0 0.00000000 Low 0
## 113 7/15/21 8/19/21 3 0 45 1.50000000 Low 48
## 114 7/15/21 8/19/21 0 5 46 1.53333333 Low 51
## 115 7/15/21 8/19/21 0 0 13 0.43333333 Low 13
## 116 7/15/21 8/19/21 0 1 60 2.00000000 Low 61
## 117 7/15/21 8/19/21 0 1 32 1.06666667 Low 33
## 118 7/15/21 8/19/21 0 2 71 2.36666667 Low 73
## 119 7/15/21 8/19/21 0 1 68 2.26666667 Low 69
## 120 7/8/21 8/12/21 0 1 51 1.70000000 Low 52
## 121 7/8/21 8/12/21 1 0 102 3.40000000 Low 103
## 122 7/8/21 8/12/21 0 0 16 0.53333333 Low 16
## 123 7/8/21 8/12/21 0 0 42 1.40000000 Low 42
## 124 7/15/21 8/19/21 0 0 59 1.96666667 Low 59
## 125 7/15/21 8/19/21 0 0 18 0.60000000 Low 18
## 126 7/15/21 8/19/21 2 0 98 3.26666667 Low 100
## 127 7/15/21 8/19/21 0 0 54 1.80000000 Low 54
## 128 7/15/21 8/19/21 5 1 102 3.40000000 Low 108
## 129 7/15/21 8/19/21 0 2 62 2.06666667 Low 64
## 130 7/22/21 8/26/21 1 3 92 3.06666667 Low 96
## 131 7/22/21 8/26/21 0 0 5 0.16666667 Low 5
## 132 7/8/21 8/12/21 0 2 72 2.40000000 Low 74
## 133 7/8/21 8/12/21 1 1 76 2.53333333 Low 78
## 134 7/8/21 8/12/21 1 0 40 1.33333333 Low 41
## 135 7/22/21 8/26/21 0 0 0 0.00000000 Low 0
## 136 7/22/21 8/26/21 2 0 39 1.30000000 Low 41
## 137 7/22/21 8/26/21 0 1 42 1.40000000 Low 43
## 138 7/22/21 8/26/21 0 0 38 1.26666667 Low 38
## 139 7/22/21 8/26/21 2 0 34 1.13333333 Low 36
## 140 7/22/21 8/26/21 4 1 75 2.50000000 Low 80
## 141 7/8/21 8/12/21 1 0 48 1.60000000 Low 49
## 142 7/8/21 8/12/21 2 2 40 1.33333333 Low 44
## 143 7/8/21 8/12/21 3 4 55 1.83333333 Low 62
## 144 7/8/21 8/12/21 1 5 36 1.20000000 Low 42
## 145 7/8/21 8/12/21 1 0 14 0.46666667 Low 15
## 146 7/8/21 8/12/21 0 0 55 1.83333333 Low 55
## 147 7/8/21 8/12/21 0 0 28 0.93333333 Low 28
## 148 7/8/21 8/12/21 1 2 55 1.83333333 Low 58
## 149 7/8/21 8/12/21 0 0 43 1.43333333 Low 43
## 150 7/8/21 8/12/21 2 1 45 1.50000000 Low 48
## 151 7/8/21 8/12/21 39 6 32 1.06666667 Low 77
## 152 7/8/21 8/12/21 1 2 65 2.16666667 Low 68
## 153 7/8/21 8/12/21 1 0 50 1.66666667 Low 51
## 154 7/8/21 8/12/21 1 1 28 0.93333333 Low 30
## 155 7/8/21 8/12/21 1 0 51 1.70000000 Low 52
## 156 7/8/21 8/12/21 1 0 55 1.83333333 Low 56
## 157 7/8/21 8/12/21 6 1 32 1.06666667 Medium 39
## 158 7/15/21 8/19/21 0 2 43 1.43333333 Medium 45
## 159 7/15/21 8/19/21 0 0 63 2.10000000 Medium 63
## 160 7/15/21 8/19/21 3 0 78 2.60000000 Medium 81
## 161 7/15/21 8/19/21 0 0 60 2.00000000 Medium 60
## 162 7/15/21 8/19/21 0 0 61 2.03333333 Medium 61
## 163 7/15/21 8/19/21 0 0 35 1.16666667 Medium 35
## 164 7/15/21 8/19/21 0 1 106 3.53333333 Medium 107
## 165 7/15/21 8/19/21 3 1 83 2.76666667 Medium 87
## 166 7/15/21 8/19/21 1 1 92 3.06666667 Medium 94
## 167 7/15/21 8/19/21 1 1 26 0.86666667 Medium 28
## 168 7/15/21 8/19/21 1 1 54 1.80000000 Medium 56
## 169 7/15/21 8/19/21 0 0 20 0.66666667 Medium 20
## 170 7/15/21 8/19/21 0 0 41 1.36666667 Medium 41
## 171 7/15/21 8/19/21 1 0 16 0.53333333 Medium 17
## 172 7/22/21 8/26/21 0 0 9 0.30000000 Medium 9
## 173 7/22/21 8/26/21 0 0 0 0.00000000 Medium 0
## 174 7/22/21 8/26/21 0 0 9 0.30000000 Medium 9
## 175 7/22/21 8/26/21 0 0 0 0.00000000 Medium 0
## 176 7/22/21 8/26/21 0 1 5 0.16666667 Medium 6
## 177 7/8/21 8/12/21 0 3 77 2.56666667 Medium 80
## 178 7/8/21 8/12/21 0 0 76 2.53333333 Medium 76
## 179 7/8/21 8/12/21 0 0 30 1.00000000 Medium 30
## 180 7/22/21 8/26/21 18 4 66 2.20000000 Medium 88
## 181 7/22/21 8/26/21 2 0 58 1.93333333 Medium 60
## 182 7/22/21 8/26/21 2 0 98 3.26666667 Medium 100
## 183 7/22/21 8/26/21 0 0 37 1.23333333 Medium 37
## 184 7/8/21 8/12/21 6 0 80 2.66666667 Medium 86
## 185 7/8/21 8/12/21 0 0 42 1.40000000 Medium 42
## 186 7/8/21 8/12/21 0 5 31 1.03333333 Medium 36
## 187 7/8/21 8/12/21 0 2 36 1.20000000 Medium 38
## 188 7/8/21 8/12/21 0 0 1 0.03333333 Medium 1
## 189 7/8/21 8/12/21 0 0 23 0.76666667 Medium 23
## 190 7/8/21 8/12/21 0 2 64 2.13333333 Medium 66
## 191 7/8/21 8/12/21 1 2 32 1.06666667 Medium 35
## 192 7/8/21 8/12/21 1 2 47 1.56666667 Medium 50
## 193 7/8/21 8/12/21 1 2 121 4.03333333 Medium 124
## 194 7/8/21 8/12/21 0 3 79 2.63333333 Medium 82
## 195 7/8/21 8/12/21 1 0 42 1.40000000 Medium 43
## 196 7/8/21 8/12/21 0 2 68 2.26666667 Medium 70
## 197 7/8/21 8/12/21 2 1 8 0.26666667 Medium 11
## 198 7/8/21 8/12/21 1 0 44 1.46666667 Medium 45
## 199 7/15/21 8/19/21 0 2 55 1.83333333 Medium 57
## 200 7/15/21 8/19/21 0 2 63 2.10000000 Medium 65
## 201 7/15/21 8/19/21 0 1 45 1.50000000 Medium 46
## 202 7/15/21 8/19/21 0 2 48 1.60000000 Medium 50
## 203 <NA> <NA> NA NA NA NA <NA> NA
## 204 <NA> <NA> NA NA NA NA <NA> NA
## 205 <NA> <NA> NA NA NA NA <NA> NA
## 206 <NA> <NA> NA NA NA NA <NA> NA
## 207 <NA> <NA> NA NA NA NA <NA> NA
## 208 <NA> <NA> NA NA NA NA <NA> NA
## 209 <NA> <NA> NA NA NA NA <NA> NA
## 210 <NA> <NA> NA NA NA NA <NA> NA
## 211 <NA> <NA> NA NA NA NA <NA> NA
## 212 <NA> <NA> NA NA NA NA <NA> NA
## 213 <NA> <NA> NA NA NA NA <NA> NA
## p_adults p_pupae p_larvae He_remain He_remain_exp He_remain_end
## 1 0.8899083 0.082568807 0.027522936 0.9986008 0.9719424 0.9705825
## 2 0.8000000 0.100000000 0.100000000 0.9986008 0.9719424 0.9705825
## 3 0.9719626 0.018691589 0.009345794 0.9986008 0.9825008 0.9811261
## 4 0.9655172 0.011494253 0.022988506 0.9986008 0.9825008 0.9811261
## 5 0.9767442 0.000000000 0.023255814 0.9986008 0.9825008 0.9811261
## 6 0.9880952 0.000000000 0.011904762 0.9986008 0.9806907 0.9793186
## 7 0.9600000 0.000000000 0.040000000 0.9986008 0.9806907 0.9793186
## 8 1.0000000 0.000000000 0.000000000 0.9986008 0.9806907 0.9793186
## 9 0.9772727 0.000000000 0.022727273 0.9986008 0.9806907 0.9793186
## 10 0.9545455 0.000000000 0.045454545 0.9986008 0.9386647 0.9373514
## 11 0.8769231 0.000000000 0.123076923 0.9986008 0.9386647 0.9373514
## 12 0.9821429 0.008928571 0.008928571 0.9986008 0.9797579 0.9783870
## 13 0.9782609 0.000000000 0.021739130 0.9986008 0.9797579 0.9783870
## 14 1.0000000 0.000000000 0.000000000 0.9986008 0.9797579 0.9783870
## 15 0.9473684 0.000000000 0.052631579 0.9986008 0.9797579 0.9783870
## 16 0.9491525 0.033898305 0.016949153 0.9986008 0.9797579 0.9783870
## 17 0.9821429 0.017857143 0.000000000 0.9986008 0.9513318 0.9500007
## 18 0.9090909 0.036363636 0.054545455 0.9986008 0.9803784 0.9790066
## 19 1.0000000 0.000000000 0.000000000 0.9986008 0.9803784 0.9790066
## 20 0.9375000 0.062500000 0.000000000 0.9986008 0.9803784 0.9790066
## 21 0.8800000 0.080000000 0.040000000 0.9986008 0.9836255 0.9822492
## 22 0.9893617 0.010638298 0.000000000 0.9986008 0.9836255 0.9822492
## 23 0.9824561 0.008771930 0.008771930 0.9986008 0.9836255 0.9822492
## 24 0.9763314 0.000000000 0.023668639 0.9986008 0.9802809 0.9789093
## 25 0.9864865 0.013513514 0.000000000 0.9986008 0.9802809 0.9789093
## 26 0.9390244 0.000000000 0.060975610 0.9986008 0.9827964 0.9814213
## 27 0.9767442 0.011627907 0.011627907 0.9986008 0.9827964 0.9814213
## 28 1.0000000 0.000000000 0.000000000 0.9986008 0.9827964 0.9814213
## 29 0.9411765 0.019607843 0.039215686 0.9986008 0.9827964 0.9814213
## 30 1.0000000 0.000000000 0.000000000 0.9986008 0.9827964 0.9814213
## 31 1.0000000 0.000000000 0.000000000 0.9986008 0.9747167 0.9733529
## 32 0.9450549 0.043956044 0.010989011 0.9986008 0.9747167 0.9733529
## 33 0.9680000 0.016000000 0.016000000 0.9986008 0.9747167 0.9733529
## 34 0.9574468 0.028368794 0.014184397 0.9986008 0.9739310 0.9725683
## 35 1.0000000 0.000000000 0.000000000 0.9986008 0.9739310 0.9725683
## 36 0.9838710 0.000000000 0.016129032 0.9986008 0.9852490 0.9838705
## 37 1.0000000 0.000000000 0.000000000 0.9986008 0.9852490 0.9838705
## 38 1.0000000 0.000000000 0.000000000 0.9986008 0.9852490 0.9838705
## 39 1.0000000 0.000000000 0.000000000 0.9986008 0.9852490 0.9838705
## 40 0.9263158 0.052631579 0.021052632 0.9986008 0.9822146 0.9808403
## 41 0.8750000 0.025000000 0.100000000 0.9986008 0.9822146 0.9808403
## 42 0.9702970 0.009900990 0.019801980 0.9986008 0.9855751 0.9841961
## 43 0.9220779 0.064935065 0.012987013 0.9986008 0.9855751 0.9841961
## 44 0.9516129 0.016129032 0.032258065 0.9986008 0.9855751 0.9841961
## 45 0.9347826 0.021739130 0.043478261 0.9986008 0.9855751 0.9841961
## 46 0.9891304 0.000000000 0.010869565 0.9986008 0.9855751 0.9841961
## 47 1.0000000 0.000000000 0.000000000 0.9986008 0.9287849 0.9274854
## 48 0.9102564 0.038461538 0.051282051 0.9986008 0.9751338 0.9737695
## 49 0.6403509 0.035087719 0.324561404 0.9986008 0.9751338 0.9737695
## 50 0.9605263 0.026315789 0.013157895 0.9986008 0.9751338 0.9737695
## 51 0.9021739 0.076086957 0.021739130 0.9986008 0.9751338 0.9737695
## 52 1.0000000 0.000000000 0.000000000 0.9986008 0.9829430 0.9815677
## 53 0.9809524 0.009523810 0.009523810 0.9986008 0.9829430 0.9815677
## 54 0.9420290 0.043478261 0.014492754 0.9986008 0.9829430 0.9815677
## 55 0.9868421 0.000000000 0.013157895 0.9986008 0.9829430 0.9815677
## 56 0.9923077 0.000000000 0.007692308 0.9986008 0.9829430 0.9815677
## 57 1.0000000 0.000000000 0.000000000 0.9986008 0.9829430 0.9815677
## 58 0.9920000 0.000000000 0.008000000 0.9986008 0.9794514 0.9780809
## 59 0.8181818 0.000000000 0.181818182 0.9986008 0.9794514 0.9780809
## 60 1.0000000 0.000000000 0.000000000 0.9986008 0.9794514 0.9780809
## 61 0.9645390 0.014184397 0.021276596 0.9986008 0.9158403 0.9145589
## 62 0.9862069 0.000000000 0.013793103 0.9986008 0.9158403 0.9145589
## 63 0.9531250 0.031250000 0.015625000 0.9986008 0.9428291 0.9415099
## 64 0.9473684 0.000000000 0.052631579 0.9986008 0.9428291 0.9415099
## 65 1.0000000 0.000000000 0.000000000 0.9986008 0.9718326 0.9704729
## 66 1.0000000 0.000000000 0.000000000 0.9986008 0.9718326 0.9704729
## 67 0.9859155 0.014084507 0.000000000 0.9986008 0.9718326 0.9704729
## 68 0.9805825 0.009708738 0.009708738 0.9986008 0.9718326 0.9704729
## 69 1.0000000 0.000000000 0.000000000 0.9986008 0.9718326 0.9704729
## 70 0.9901961 0.000000000 0.009803922 0.9986008 0.9809878 0.9796153
## 71 1.0000000 0.000000000 0.000000000 0.9986008 0.9809878 0.9796153
## 72 1.0000000 0.000000000 0.000000000 0.9986008 0.9809878 0.9796153
## 73 0.9868421 0.000000000 0.013157895 0.9986008 0.9809878 0.9796153
## 74 1.0000000 0.000000000 0.000000000 0.9986008 0.9809878 0.9796153
## 75 1.0000000 0.000000000 0.000000000 0.9986008 0.9809878 0.9796153
## 76 1.0000000 0.000000000 0.000000000 0.9986008 0.9835482 0.9821721
## 77 1.0000000 0.000000000 0.000000000 0.9986008 0.9835482 0.9821721
## 78 1.0000000 0.000000000 0.000000000 0.9986008 0.9835482 0.9821721
## 79 1.0000000 0.000000000 0.000000000 0.9986008 0.9835482 0.9821721
## 80 1.0000000 0.000000000 0.000000000 0.9986008 0.9835482 0.9821721
## 81 0.9777778 0.022222222 0.000000000 0.9986008 0.9855894 0.9842104
## 82 1.0000000 0.000000000 0.000000000 0.9986008 0.9855894 0.9842104
## 83 1.0000000 0.000000000 0.000000000 0.9986008 0.9855894 0.9842104
## 84 0.9878049 0.000000000 0.012195122 0.9986008 0.9855894 0.9842104
## 85 1.0000000 0.000000000 0.000000000 0.9986008 0.9855894 0.9842104
## 86 1.0000000 0.000000000 0.000000000 0.9986008 0.9855894 0.9842104
## 87 0.9729730 0.027027027 0.000000000 0.9986008 0.9855894 0.9842104
## 88 0.9838710 0.000000000 0.016129032 0.9986008 0.9717830 0.9704233
## 89 0.9444444 0.000000000 0.055555556 0.9986008 0.9717830 0.9704233
## 90 0.9565217 0.014492754 0.028985507 0.9986008 0.9717830 0.9704233
## 91 0.9901961 0.009803922 0.000000000 0.9986008 0.9682092 0.9668545
## 92 0.9509804 0.019607843 0.029411765 0.9986008 0.9682092 0.9668545
## 93 0.9157895 0.084210526 0.000000000 0.9986008 0.9682092 0.9668545
## 94 1.0000000 0.000000000 0.000000000 0.9986008 0.9682092 0.9668545
## 95 0.9200000 0.080000000 0.000000000 0.9986008 0.9682092 0.9668545
## 96 1.0000000 0.000000000 0.000000000 0.9986008 0.9682092 0.9668545
## 97 1.0000000 0.000000000 0.000000000 0.5456041 0.9801761 0.5347881
## 98 0.9777778 0.000000000 0.022222222 0.5456041 0.9801761 0.5347881
## 99 1.0000000 0.000000000 0.000000000 0.5456041 0.9801761 0.5347881
## 100 1.0000000 0.000000000 0.000000000 0.5456041 0.9801761 0.5347881
## 101 0.9661017 0.016949153 0.016949153 0.5456041 0.9801761 0.5347881
## 102 0.9655172 0.034482759 0.000000000 0.5411689 0.9684567 0.5240986
## 103 1.0000000 0.000000000 0.000000000 0.5375090 0.9279284 0.4987699
## 104 0.9827586 0.000000000 0.017241379 0.5487325 0.9725983 0.5336963
## 105 0.9750000 0.008333333 0.016666667 0.5487325 0.9725983 0.5336963
## 106 0.9166667 0.011904762 0.071428571 0.5546865 0.9630618 0.5341974
## 107 0.9896907 0.010309278 0.000000000 0.5507359 0.9403195 0.5178677
## 108 0.9578947 0.021052632 0.021052632 0.5507359 0.9403195 0.5178677
## 109 0.9714286 0.009523810 0.019047619 0.5507359 0.9403195 0.5178677
## 110 0.9473684 0.039473684 0.013157895 0.5532470 0.9656184 0.5342255
## 111 0.9272727 0.018181818 0.054545455 0.5532470 0.9656184 0.5342255
## 112 NaN NaN NaN 0.5544299 0.6449276 0.3575672
## 113 0.9375000 0.000000000 0.062500000 0.5071881 0.9790690 0.4965722
## 114 0.9019608 0.098039216 0.000000000 0.5071881 0.9790690 0.4965722
## 115 1.0000000 0.000000000 0.000000000 0.4987816 0.9782357 0.4879260
## 116 0.9836066 0.016393443 0.000000000 0.4987816 0.9782357 0.4879260
## 117 0.9696970 0.030303030 0.000000000 0.4987816 0.9782357 0.4879260
## 118 0.9726027 0.027397260 0.000000000 0.4987816 0.9782357 0.4879260
## 119 0.9855072 0.014492754 0.000000000 0.5541380 0.9764292 0.5410765
## 120 0.9807692 0.019230769 0.000000000 0.5448947 0.9741394 0.5308034
## 121 0.9902913 0.000000000 0.009708738 0.5448947 0.9741394 0.5308034
## 122 1.0000000 0.000000000 0.000000000 0.5448947 0.9741394 0.5308034
## 123 1.0000000 0.000000000 0.000000000 0.5448947 0.9741394 0.5308034
## 124 1.0000000 0.000000000 0.000000000 0.5284581 0.9093326 0.4805442
## 125 1.0000000 0.000000000 0.000000000 0.5529282 0.8511338 0.4706159
## 126 0.9800000 0.000000000 0.020000000 0.5435302 0.9711223 0.5278343
## 127 1.0000000 0.000000000 0.000000000 0.5476417 0.9565572 0.5238506
## 128 0.9444444 0.009259259 0.046296296 0.5430646 0.9758713 0.5299612
## 129 0.9687500 0.031250000 0.000000000 0.5430646 0.9758713 0.5299612
## 130 0.9583333 0.031250000 0.010416667 0.5509294 0.9167235 0.5050499
## 131 1.0000000 0.000000000 0.000000000 0.5431402 0.7694643 0.4179270
## 132 0.9729730 0.027027027 0.000000000 0.5460067 0.9556060 0.5217673
## 133 0.9743590 0.012820513 0.012820513 0.5460067 0.9556060 0.5217673
## 134 0.9756098 0.000000000 0.024390244 0.5460067 0.9556060 0.5217673
## 135 NaN NaN NaN 0.5523333 0.6083089 0.3359893
## 136 0.9512195 0.000000000 0.048780488 0.5463114 0.9631949 0.5262044
## 137 0.9767442 0.023255814 0.000000000 0.5463114 0.9631949 0.5262044
## 138 1.0000000 0.000000000 0.000000000 0.5463114 0.9631949 0.5262044
## 139 0.9444444 0.000000000 0.055555556 0.5463114 0.9631949 0.5262044
## 140 0.9375000 0.012500000 0.050000000 0.5542119 0.8901046 0.4933065
## 141 0.9795918 0.000000000 0.020408163 0.5489205 0.9707940 0.5328887
## 142 0.9090909 0.045454545 0.045454545 0.5489205 0.9707940 0.5328887
## 143 0.8870968 0.064516129 0.048387097 0.5489205 0.9707940 0.5328887
## 144 0.8571429 0.119047619 0.023809524 0.5489205 0.9707940 0.5328887
## 145 0.9333333 0.000000000 0.066666667 0.5489205 0.9707940 0.5328887
## 146 1.0000000 0.000000000 0.000000000 0.5419774 0.9678511 0.5245534
## 147 1.0000000 0.000000000 0.000000000 0.5419774 0.9678511 0.5245534
## 148 0.9482759 0.034482759 0.017241379 0.5419774 0.9678511 0.5245534
## 149 1.0000000 0.000000000 0.000000000 0.5541201 0.9734332 0.5393989
## 150 0.9375000 0.020833333 0.041666667 0.5541201 0.9734332 0.5393989
## 151 0.4155844 0.077922078 0.506493506 0.5541201 0.9734332 0.5393989
## 152 0.9558824 0.029411765 0.014705882 0.5541201 0.9734332 0.5393989
## 153 0.9803922 0.000000000 0.019607843 0.5274650 0.9808799 0.5173798
## 154 0.9333333 0.033333333 0.033333333 0.5274650 0.9808799 0.5173798
## 155 0.9807692 0.000000000 0.019230769 0.5274650 0.9808799 0.5173798
## 156 0.9821429 0.000000000 0.017857143 0.5274650 0.9808799 0.5173798
## 157 0.8205128 0.025641026 0.153846154 0.6596334 0.9422801 0.6215595
## 158 0.9555556 0.044444444 0.000000000 0.6754272 0.9327651 0.6300149
## 159 1.0000000 0.000000000 0.000000000 0.6803810 0.9756391 0.6638064
## 160 0.9629630 0.000000000 0.037037037 0.6803810 0.9756391 0.6638064
## 161 1.0000000 0.000000000 0.000000000 0.6803810 0.9756391 0.6638064
## 162 1.0000000 0.000000000 0.000000000 0.6803810 0.9756391 0.6638064
## 163 1.0000000 0.000000000 0.000000000 0.6831935 0.8504137 0.5809971
## 164 0.9906542 0.009345794 0.000000000 0.6807728 0.9721768 0.6618315
## 165 0.9540230 0.011494253 0.034482759 0.6807728 0.9721768 0.6618315
## 166 0.9787234 0.010638298 0.010638298 0.6807728 0.9721768 0.6618315
## 167 0.9285714 0.035714286 0.035714286 0.6620374 0.9688824 0.6414364
## 168 0.9642857 0.017857143 0.017857143 0.6620374 0.9688824 0.6414364
## 169 1.0000000 0.000000000 0.000000000 0.6620374 0.9688824 0.6414364
## 170 1.0000000 0.000000000 0.000000000 0.6620374 0.9688824 0.6414364
## 171 0.9411765 0.000000000 0.058823529 0.6620374 0.9688824 0.6414364
## 172 1.0000000 0.000000000 0.000000000 0.6827942 0.9861527 0.6733394
## 173 NaN NaN NaN 0.6827942 0.9861527 0.6733394
## 174 1.0000000 0.000000000 0.000000000 0.6827942 0.9861527 0.6733394
## 175 NaN NaN NaN 0.6827942 0.9861527 0.6733394
## 176 0.8333333 0.166666667 0.000000000 0.6827298 0.8930853 0.6097360
## 177 0.9625000 0.037500000 0.000000000 0.6799082 0.9736031 0.6619607
## 178 1.0000000 0.000000000 0.000000000 0.6799082 0.9736031 0.6619607
## 179 1.0000000 0.000000000 0.000000000 0.6799082 0.9736031 0.6619607
## 180 0.7500000 0.045454545 0.204545455 0.6804075 0.9418041 0.6408105
## 181 0.9666667 0.000000000 0.033333333 0.6768800 0.9677851 0.6550744
## 182 0.9800000 0.000000000 0.020000000 0.6802035 0.9422407 0.6409154
## 183 1.0000000 0.000000000 0.000000000 0.6802035 0.9422407 0.6409154
## 184 0.9302326 0.000000000 0.069767442 0.6819561 0.9777829 0.6668050
## 185 1.0000000 0.000000000 0.000000000 0.6819561 0.9777829 0.6668050
## 186 0.8611111 0.138888889 0.000000000 0.6819561 0.9777829 0.6668050
## 187 0.9473684 0.052631579 0.000000000 0.6826989 0.9622026 0.6568947
## 188 1.0000000 0.000000000 0.000000000 0.6807065 0.4714889 0.3209456
## 189 1.0000000 0.000000000 0.000000000 0.6819044 0.9771934 0.6663525
## 190 0.9696970 0.030303030 0.000000000 0.6819044 0.9771934 0.6663525
## 191 0.9142857 0.057142857 0.028571429 0.6819044 0.9771934 0.6663525
## 192 0.9400000 0.040000000 0.020000000 0.6800863 0.9695908 0.6594054
## 193 0.9758065 0.016129032 0.008064516 0.6800863 0.9695908 0.6594054
## 194 0.9634146 0.036585366 0.000000000 0.6820356 0.9822570 0.6699342
## 195 0.9767442 0.000000000 0.023255814 0.6820356 0.9822570 0.6699342
## 196 0.9714286 0.028571429 0.000000000 0.6820356 0.9822570 0.6699342
## 197 0.7272727 0.090909091 0.181818182 0.6820356 0.9822570 0.6699342
## 198 0.9777778 0.000000000 0.022222222 0.6820356 0.9822570 0.6699342
## 199 0.9649123 0.035087719 0.000000000 0.6804431 0.9819642 0.6681708
## 200 0.9692308 0.030769231 0.000000000 0.6804431 0.9819642 0.6681708
## 201 0.9782609 0.021739130 0.000000000 0.6804431 0.9819642 0.6681708
## 202 0.9600000 0.040000000 0.000000000 0.6804431 0.9819642 0.6681708
## 203 NA NA NA 0.5367998 0.8056432 0.4324691
## 204 NA NA NA 0.5386585 0.7452523 0.4014365
## 205 NA NA NA 0.5540444 0.4950540 0.2742819
## 206 NA NA NA 0.5532775 0.7440450 0.4116634
## 207 NA NA NA 0.5528880 0.9892872 0.5469651
## 208 NA NA NA 0.5427371 0.4855518 0.2635269
## 209 NA NA NA 0.5518545 0.9696290 0.5350941
## 210 NA NA NA 0.6802331 0.9285014 0.6315974
## 211 NA NA NA 0.6825509 0.8944237 0.6104897
## 212 NA NA NA 0.6819650 0.4841994 0.3302071
## 213 NA NA NA 0.6809168 0.9614965 0.6546991SUM_POP_He_ALL_extinct<- Rmisc::summarySE(data_phenotyping_withextinct,
measurevar = c("Lambda"),
groupvars = c("Genetic_Diversity","He_remain", "He_remain_end",
"Population"),
na.rm=TRUE)
#Add extinction variable
SUM_POP_He_ALL_extinct$Extinction_finalstatus <-ifelse(is.na(SUM_POP_He_ALL_extinct$Lambda),
"yes","no")
#Replace extinct by NA
SUM_POP_He_ALL_extinct$Lambda <- ifelse(is.na(SUM_POP_He_ALL_extinct$Lambda),
0,SUM_POP_He_ALL_extinct$Lambda)
SUM_POP_He_ALL_extinct$Genetic_Diversity <- ifelse(is.na(SUM_POP_He_ALL_extinct$Genetic_Diversity),
SUM_POP_He_ALL_extinct$Population,
SUM_POP_He_ALL_extinct$Genetic_Diversity)
SUM_POP_He_ALL_extinct$Genetic_Diversity[grep("Low",as.character(SUM_POP_He_ALL_extinct$Population))] <- "Low"
SUM_POP_He_ALL_extinct$Genetic_Diversity[grep("Med",as.character(SUM_POP_He_ALL_extinct$Population))] <- "Medium"
SUM_POP_He_ALL_extinct$Genetic_Diversity[grep("High",as.character(SUM_POP_He_ALL_extinct$Population))] <- "High"
#Order levels
SUM_POP_He_ALL_extinct$Genetic_Diversity <- factor(SUM_POP_He_ALL_extinct$Genetic_Diversity, levels = c("High", "Medium", "Low"))
PLOT_He_withextinct <- ggplot2::ggplot() +
geom_point(data = SUM_POP_He_ALL_extinct,
aes(x = He_remain_end, y = Lambda,
colour = Genetic_Diversity,
shape = Extinction_finalstatus,
size = N), alpha = 0.7) +
ylab("Growth rate") +
xlab("Remaining heterozygosity") +
scale_shape_manual(values = c(19, 1)) +
theme_LO +
labs(color="Genetic diversity") +
labs(shape="Extinct population") +
labs(size="Replicates")
PLOT_He_withextinct
# cowplot::save_plot(here::here("figures","Fitness_He_final.pdf"), PLOT_He_withextinct,
# base_height = 12/cm(1), base_width = 18/cm(1), dpi = 610)3 ANALYSIS
3.1 Probability of extinction
############
###### Clean dataset
############
#Prepare clean dataset
data_proba_extinction <- data_G6[,c("Block","Population","Genetic_Diversity","Extinction_finalstatus")]
data_proba_extinction$Extinction_finalstatus <- as.factor(data_proba_extinction$Extinction_finalstatus)
data_proba_extinction$y <- ifelse(data_proba_extinction$Extinction_finalstatus == "yes", 1, 0)
############
###### Analysis
############
#Analysis
mod0 <- glm(y ~ Genetic_Diversity + Block,
data = data_proba_extinction, family = binomial)
summary(mod0)##
## Call:
## glm(formula = y ~ Genetic_Diversity + Block, family = binomial,
## data = data_proba_extinction)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.26011 -0.44258 -0.30957 -0.00005 2.47474
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -21.3144 1855.7435 -0.011 0.99084
## Genetic_DiversityMedium 18.3002 1855.7433 0.010 0.99213
## Genetic_DiversityLow 18.5062 1855.7432 0.010 0.99204
## BlockBlock4 0.7402 1.2714 0.582 0.56046
## BlockBlock5 3.0006 1.1265 2.664 0.00773 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 72.388 on 83 degrees of freedom
## Residual deviance: 46.833 on 79 degrees of freedom
## AIC: 56.833
##
## Number of Fisher Scoring iterations: 18#Test genetic diversity effect
mod1 <- glm(y ~ Block ,
data = data_proba_extinction, family = binomial)
anova(mod1, mod0, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 81 58.903
## 2 79 46.833 2 12.069 0.002394 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#We keep genetic diversity
# Perform the lrtest
(A <- logLik(mod1))## 'log Lik.' -29.45146 (df=3)(B <- logLik(mod0))## 'log Lik.' -23.41669 (df=5)#Since the logLik() function gives more information than the numeric value, use the as.numeric() function to isolate the numeric value
(teststat <- -2 * (as.numeric(A)-as.numeric(B)))## [1] 12.06954#df = 5 - 3 = 2
(p.val <- pchisq(teststat, df = 2, lower.tail = FALSE))## [1] 0.002394043lmtest::lrtest(mod1, mod0)## Likelihood ratio test
##
## Model 1: y ~ Block
## Model 2: y ~ Genetic_Diversity + Block
## #Df LogLik Df Chisq Pr(>Chisq)
## 1 3 -29.451
## 2 5 -23.417 2 12.069 0.002394 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod0 <- glm(y ~ Genetic_Diversity-1 + Block ,
data = data_proba_extinction, family = binomial)
summary(mod0)##
## Call:
## glm(formula = y ~ Genetic_Diversity - 1 + Block, family = binomial,
## data = data_proba_extinction)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.26011 -0.44258 -0.30957 -0.00005 2.47474
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## Genetic_DiversityHigh -21.3144 1855.7435 -0.011 0.99084
## Genetic_DiversityMedium -3.0142 1.1293 -2.669 0.00760 **
## Genetic_DiversityLow -2.8082 1.0659 -2.635 0.00843 **
## BlockBlock4 0.7402 1.2714 0.582 0.56046
## BlockBlock5 3.0006 1.1265 2.664 0.00773 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 116.449 on 84 degrees of freedom
## Residual deviance: 46.833 on 79 degrees of freedom
## AIC: 56.833
##
## Number of Fisher Scoring iterations: 18############
###### Predicted value
############
#Extract probability of extinction
(p_high <- plogis(-21.3144))## [1] 5.536989e-10(p_med <- plogis(-3.0142))## [1] 0.04678847(p_low <- plogis(-2.8082))## [1] 0.05688267data_predict_extinct <- data.frame(Genetic_Diversity = levels(data_proba_extinction$Genetic_Diversity),
Block = "Block4")
data_predict_extinct$predict <- predict(mod0, type = "response",
re.form = NA,
newdata = data_predict_extinct)
data_predict_extinct## Genetic_Diversity Block predict
## 1 High Block4 1.160723e-09
## 2 Medium Block4 9.329490e-02
## 3 Low Block4 1.122446e-01p_high## [1] 5.536989e-10p_med## [1] 0.04678847p_low## [1] 0.05688267############
###### Posthoc
############
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High -20.07 1855.743 Inf -4451.10 4410.961
## Medium -1.77 0.650 Inf -3.32 -0.216
## Low -1.56 0.532 Inf -2.83 -0.290
##
## Results are averaged over the levels of: Block
## Results are given on the logit (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium -18.300 1855.743 Inf -0.010 0.9999
## High - Low -18.506 1855.743 Inf -0.010 0.9999
## Medium - Low -0.206 0.751 Inf -0.274 0.9594
##
## Results are averaged over the levels of: Block
## Results are given on the log odds ratio (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates3.2 Population size G6
############
###### Analysis
############
#Analysis
# mod0 <- lme4::glmer(Nb_adults ~ Genetic_Diversity + Block + (1|Population),
# data = data_G6, family = "poisson")
mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block,
data = data_G6[data_G6$Extinction_finalstatus == "no",])
mod1 <- lm(log(Nb_adults) ~ Block ,
data = data_G6[data_G6$Extinction_finalstatus == "no",])
anova(mod0, mod1) # Effect of genetic diversity ## Analysis of Variance Table
##
## Model 1: log(Nb_adults) ~ Genetic_Diversity + Block
## Model 2: log(Nb_adults) ~ Block
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 66 29.706
## 2 68 42.118 -2 -12.412 13.789 9.914e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod0 <- lm(log(Nb_adults) ~ Genetic_Diversity + Block,
data = data_G6[data_G6$Extinction_finalstatus == "no",])
summary(mod0)##
## Call:
## lm(formula = log(Nb_adults) ~ Genetic_Diversity + Block, data = data_G6[data_G6$Extinction_finalstatus ==
## "no", ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.2417 -0.2947 0.1088 0.4381 1.2624
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.11532 0.17980 28.451 < 2e-16 ***
## Genetic_DiversityMedium -0.94813 0.20540 -4.616 1.86e-05 ***
## Genetic_DiversityLow -0.80532 0.18719 -4.302 5.72e-05 ***
## BlockBlock4 -0.02077 0.18447 -0.113 0.9107
## BlockBlock5 -0.53922 0.21414 -2.518 0.0142 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6709 on 66 degrees of freedom
## Multiple R-squared: 0.3362, Adjusted R-squared: 0.2959
## F-statistic: 8.356 on 4 and 66 DF, p-value: 1.623e-05############
###### Posthoc
############
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 4.93 0.130 66 4.61 5.25
## Medium 3.98 0.159 66 3.59 4.37
## Low 4.12 0.136 66 3.79 4.46
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.948 0.205 66 4.616 0.0001
## High - Low 0.805 0.187 66 4.302 0.0002
## Medium - Low -0.143 0.207 66 -0.690 0.7704
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates3.3 Population size over time
#######################
####################### MODEL INCLUDING GENETIC DIVERSITY WITHOUT GEN 1
#######################
###################### REMOVE GEN 1 ##########################33
#If we dont consider Gen=1
PLOT_Pop_size_average
fit <- lm(log(Nb_adults+1) ~ Generation_Eggs, data = data[data$Generation_Eggs>=2,])
#second degree
fit2 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#third degree
fit3 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,3,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#fourth degree
fit4 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,4,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#fourth degree
fit5 <- lm(log(Nb_adults+1)~stats::poly(Generation_Eggs,5,raw=TRUE), data = data[data$Generation_Eggs>=2,])
#exp
fit6 <- lm(log(Nb_adults+1)~exp(Generation_Eggs), data = data[data$Generation_Eggs>=2,])
#log
fit7 <- lm(log(Nb_adults+1)~log(Generation_Eggs), data = data[data$Generation_Eggs>=2,])
fit <- glm(Nb_adults ~ Generation_Eggs,
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE),
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit3 <- glm(Nb_adults~stats::poly(Generation_Eggs,3,raw=TRUE),
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE),
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit5 <- glm(Nb_adults~stats::poly(Generation_Eggs,5,raw=TRUE),
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit6 <- glm(Nb_adults~exp(Generation_Eggs),
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit7 <- glm(Nb_adults~log(Generation_Eggs),
data = data[data$Generation_Eggs>=2,], family = "poisson")
#generate range of 50 numbers starting from 30 and ending at 160
# xx <- seq(1,6, length=100)
# plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
# log(data[data$Generation_Eggs>=2,]$Nb_adults+1),pch=19,ylim=c(0,8))
# lines(xx, predict(fit, data.frame(Generation_Eggs=xx)), col="red")
# lines(xx, predict(fit2, data.frame(Generation_Eggs=xx)), col="green")
# lines(xx, predict(fit3, data.frame(Generation_Eggs=xx)), col="blue")
# lines(xx, predict(fit4, data.frame(Generation_Eggs=xx)), col="purple")
# lines(xx, predict(fit5, data.frame(Generation_Eggs=xx)), col="yellow")
# lines(xx, predict(fit6, data.frame(Generation_Eggs=xx)), col="orange")
# lines(xx, predict(fit7, data.frame(Generation_Eggs=xx)), col="pink")
xx <- seq(1,6, length=100)
plot(data[data$Generation_Eggs>=2,]$Generation_Eggs,
data[data$Generation_Eggs>=2,]$Nb_adults,pch=19)
lines(xx, predict(fit, data.frame(Generation_Eggs=xx),type = "response"), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx),type = "response"), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx),type = "response"), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx),type = "response"), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx),type = "response"), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx),type = "response"), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx),type = "response"), col="pink")
# Check the r squared
rsq::rsq(fit,adj=TRUE)## [1] 0.4442203rsq::rsq(fit2,adj=TRUE)## [1] 0.5938535rsq::rsq(fit3,adj=TRUE)## [1] 0.5928435rsq::rsq(fit4,adj=TRUE)## [1] 0.592052rsq::rsq(fit5,adj=TRUE)## [1] 0.592052rsq::rsq(fit6,adj=TRUE)## [1] 0.121301rsq::rsq(fit7,adj=TRUE)## [1] 0.5205104# Best model is fit 2
data$gen_square <- data$Generation_Eggs^2
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE),
data = data[data$Generation_Eggs>=2,], family = "poisson")
fit2 <- glm(Nb_adults~Generation_Eggs + gen_square,
data = data[data$Generation_Eggs>=2,], family = "poisson")
# fit2 <- lm(log(Nb_adults+1)~poly(Generation_Eggs,2,raw=TRUE), data = data[data$Generation_Eggs>=2,])
# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs + gen_square, data = data[data$Generation_Eggs>=2,])
#Add the interaction with Genetic diversity
fit2 <- glm(Nb_adults~ Generation_Eggs*Genetic_Diversity +
gen_square*Genetic_Diversity,
data = data, family = "poisson")
summary(fit2)##
## Call:
## glm(formula = Nb_adults ~ Generation_Eggs * Genetic_Diversity +
## gen_square * Genetic_Diversity, family = "poisson", data = data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -15.961 -6.480 -2.662 5.254 21.115
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.070332 0.026463 191.604 < 2e-16
## Generation_Eggs 0.156559 0.017811 8.790 < 2e-16
## Genetic_DiversityMedium -0.164697 0.042046 -3.917 8.96e-05
## Genetic_DiversityLow 0.024718 0.037632 0.657 0.511290
## gen_square -0.036751 0.002566 -14.324 < 2e-16
## Generation_Eggs:Genetic_DiversityMedium 0.080423 0.029501 2.726 0.006409
## Generation_Eggs:Genetic_DiversityLow -0.096063 0.026425 -3.635 0.000278
## Genetic_DiversityMedium:gen_square -0.035166 0.004448 -7.905 2.67e-15
## Genetic_DiversityLow:gen_square -0.009935 0.003956 -2.511 0.012038
##
## (Intercept) ***
## Generation_Eggs ***
## Genetic_DiversityMedium ***
## Genetic_DiversityLow
## gen_square ***
## Generation_Eggs:Genetic_DiversityMedium **
## Generation_Eggs:Genetic_DiversityLow ***
## Genetic_DiversityMedium:gen_square ***
## Genetic_DiversityLow:gen_square *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 39168 on 486 degrees of freedom
## Residual deviance: 30719 on 478 degrees of freedom
## (17 observations deleted due to missingness)
## AIC: 33747
##
## Number of Fisher Scoring iterations: 5####### Add genetic diversity
#Test oversdispersion
mod0 <- glm(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity +
Block ,
data = data[data$Generation_Eggs>=2,], family = "poisson")
sqrt(deviance(mod0)/(nobs(mod0)-length(coef(mod0))))## [1] 6.395475sigma(mod0)## [1] 6.395475# ccl: overdispersion
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity +
Block + (1|obs),
data = data[data$Generation_Eggs>=2,], family = "poisson")
blmeco::dispersion_glmer(mod1)## [1] 1.076696sigma(mod1)## [1] 1#Convergence issue
max(abs(with(mod1@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001## [1] 0.004551091# Test effect
mod1 <- lme4::glmer(Nb_adults~Generation_Eggs*Genetic_Diversity + gen_square*Genetic_Diversity +
Block + (1|obs),
data = data[data$Generation_Eggs>=2,], family = "poisson")
mod2 <- lme4::glmer(Nb_adults~Generation_Eggs + gen_square + Genetic_Diversity +
Block + (1|obs),
data = data[data$Generation_Eggs>=2,], family = "poisson")
anova(mod1, mod2, test ="Chisq") ## Data: data[data$Generation_Eggs >= 2, ]
## Models:
## mod2: Nb_adults ~ Generation_Eggs + gen_square + Genetic_Diversity +
## mod2: Block + (1 | obs)
## mod1: Nb_adults ~ Generation_Eggs * Genetic_Diversity + gen_square *
## mod1: Genetic_Diversity + Block + (1 | obs)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod2 8 4669.1 4701.1 -2326.6 4653.1
## mod1 12 4660.8 4708.8 -2318.4 4636.8 16.376 4 0.002554 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100),
Generation_Eggs = rep(seq(2,6, length = 100),each=3),
Block = "Block4",
gen_square = (rep(seq(2,6, length = 100),each=3))^2,
Estimates = NA)
filldata$Estimates <- predict(mod1, newdata=filldata, type = "response", re.form=~0)
#Plot
PLOT_Pop_size <- ggplot2::ggplot(data = data[data$Generation_Eggs>=2,]) +
geom_line(data = filldata, aes(x = Generation_Eggs,
y = Estimates,
colour = Genetic_Diversity), size = 1.3) +
geom_point(data = data,
aes(x = Generation_Eggs,
y = Nb_adults,
colour = Genetic_Diversity),
position = position_dodge(0.4), size = 1.4, alpha = 0.6) +
ylab("Population size") +
labs(color="Genetic diversity") +
xlab("Generation") +
theme_LO
PLOT_Pop_size
emmeans::emmeans(mod1, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 4.88 0.0790 Inf 4.69 5.07
## Medium 4.19 0.0885 Inf 3.98 4.40
## Low 4.04 0.0733 Inf 3.87 4.22
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.688 0.119 Inf 5.786 <.0001
## High - Low 0.835 0.107 Inf 7.778 <.0001
## Medium - Low 0.147 0.115 Inf 1.287 0.4028
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates###################################################################
###################################################################
###################################################################
#################### INCLUDING GEN 1 ###########################
PLOT_Pop_size_average
fit <- glm(Nb_adults ~ Generation_Eggs, data = data, family = "poisson")
#second degree
fit2 <- glm(Nb_adults~stats::poly(Generation_Eggs,2,raw=TRUE), data = data, family = "poisson")
#third degree
fit3 <- glm(Nb_adults~stats::poly(Generation_Eggs,3,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")
#fourth degree
fit5 <- glm(Nb_adults~stats::poly(Generation_Eggs,5,raw=TRUE), data = data, family = "poisson")
#exp
fit6 <- glm(Nb_adults~exp(Generation_Eggs), data = data, family = "poisson")
#log
fit7 <- glm(Nb_adults~log(Generation_Eggs), data = data, family = "poisson")
#generate range of 50 numbers starting from 30 and ending at 160
xx <- seq(1,6, length=100)
plot(data$Generation_Eggs,
data$Nb_adults,pch=19)
lines(xx, predict(fit, data.frame(Generation_Eggs=xx),type = "response"), col="red")
lines(xx, predict(fit2, data.frame(Generation_Eggs=xx),type = "response"), col="green")
lines(xx, predict(fit3, data.frame(Generation_Eggs=xx),type = "response"), col="blue")
lines(xx, predict(fit4, data.frame(Generation_Eggs=xx),type = "response"), col="purple")
lines(xx, predict(fit5, data.frame(Generation_Eggs=xx),type = "response"), col="yellow")
lines(xx, predict(fit6, data.frame(Generation_Eggs=xx),type = "response"), col="orange")
lines(xx, predict(fit7, data.frame(Generation_Eggs=xx),type = "response"), col="pink")
# Check the r squared
summary(fit)$adj.r.squared## NULLsummary(fit2)$adj.r.squared## NULLsummary(fit3)$adj.r.squared## NULLsummary(fit4)$adj.r.squared # best r square## NULLsummary(fit5)$adj.r.squared ## NULLsummary(fit6)$adj.r.squared ## NULLsummary(fit7)$adj.r.squared ## NULLrsq::rsq(fit,adj=TRUE)## [1] 0.1067643rsq::rsq(fit2,adj=TRUE)## [1] 0.1486733rsq::rsq(fit3,adj=TRUE)## [1] 0.5240033rsq::rsq(fit4,adj=TRUE)## [1] 0.5980368rsq::rsq(fit5,adj=TRUE)## [1] 0.5989426rsq::rsq(fit6,adj=TRUE)## [1] 0.07185974rsq::rsq(fit7,adj=TRUE)## [1] 0.04937507#
fit4 <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE), data = data, family = "poisson")
#######################
####################### MODEL INCLUDING GENETIC DIVERSITY and GEN 1
#######################
#Test oversdispersion
fit4 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity + (1|Block), data = data, family = "poisson")
blmeco::dispersion_glmer(fit4)## [1] 5.760264mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity +
(1|Block) + (1|obs), data = data, family = "poisson")
blmeco::dispersion_glmer(mod)## [1] 1.083377#Convergence issue
max(abs(with(mod@optinfo$derivs, solve(Hessian, gradient)))) #should be <0.001## [1] 0.06384154#Predict data
filldata <- data.frame(Genetic_Diversity = rep(levels(data$Genetic_Diversity),100),
Generation_Eggs = rep(seq(1,6, length = 100),each=3),
Block = "Block4",
Estimates = NA)
filldata$Estimates <- predict(fit4, newdata=filldata, type = "response")
#Plot
PLOT_Pop_size <- ggplot2::ggplot(data = data) +
geom_line(data = filldata, aes(x = Generation_Eggs,
y = Estimates,
colour = Genetic_Diversity), size = 1.3) +
geom_point(data = data,
aes(x = Generation_Eggs,
y = Nb_adults,
colour = Genetic_Diversity),
position = position_dodge(0.4), size = 1.4, alpha = 0.6) +
ylab("Population size") +
labs(color="Genetic diversity") +
xlab("Generation") +
theme_LO
PLOT_Pop_size
# fit2 <- lm(log(Nb_adults+1)~ Generation_Eggs*Genetic_Diversity +
# gen_square*Genetic_Diversity,
# data = data[data$Generation_Eggs>=2,])
# fit2_test <- lm(log(Nb_adults+1)~ Generation_Eggs+Genetic_Diversity +
# gen_square,
# data = data[data$Generation_Eggs>=2,])
mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity +
(1|Block) + (1|obs), data = data, family = "poisson")
mod1 <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity +
(1|Block) + (1|obs), data = data, family = "poisson")
anova(mod,mod1, test = "Chisq")## Data: data
## Models:
## mod1: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity +
## mod1: (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +
## mod: (1 | Block) + (1 | obs)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 9 5611.2 5648.9 -2796.6 5593.2
## mod 17 5594.8 5666.0 -2780.4 5560.8 32.47 8 7.672e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#Signifcant interaction
summary(mod)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +
## (1 | Block) + (1 | obs)
## Data: data
##
## AIC BIC logLik deviance df.resid
## 5594.8 5666.0 -2780.4 5560.8 470
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.03685 -0.03969 0.00800 0.05695 0.26801
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.64646 0.8040
## Block (Intercept) 0.07945 0.2819
## Number of obs: 487, groups: obs, 487; Block, 3
##
## Fixed effects:
## Estimate
## (Intercept) -1.993477
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 10.858800
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -5.147603
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 0.942314
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -0.058998
## Genetic_DiversityMedium -0.352704
## Genetic_DiversityLow -0.787854
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.772686
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium -0.503761
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.095237
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium -0.006018
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 1.489498
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow -0.833931
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.138807
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow -0.007028
## Std. Error
## (Intercept) 1.238244
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 1.985855
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 1.032426
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 0.213266
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 0.015171
## Genetic_DiversityMedium 1.833248
## Genetic_DiversityLow 1.630486
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 2.976126
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 1.553127
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.322073
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.022992
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 2.637956
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 1.373155
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.284239
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.020268
## z value
## (Intercept) -1.610
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 5.468
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -4.986
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 4.418
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -3.889
## Genetic_DiversityMedium -0.192
## Genetic_DiversityLow -0.483
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.260
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium -0.324
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.296
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium -0.262
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 0.565
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow -0.607
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.488
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow -0.347
## Pr(>|z|)
## (Intercept) 0.107415
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 4.55e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 6.17e-07
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 9.94e-06
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 0.000101
## Genetic_DiversityMedium 0.847434
## Genetic_DiversityLow 0.628952
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.795151
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 0.745671
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.767459
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.793538
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 0.572318
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 0.543645
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.625306
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.728786
##
## (Intercept)
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 ***
## Genetic_DiversityMedium
## Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 0.398806 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?summary(mod1)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity +
## (1 | Block) + (1 | obs)
## Data: data
##
## AIC BIC logLik deviance df.resid
## 5611.2 5648.9 -2796.6 5593.2 478
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.05397 -0.05185 0.01426 0.07054 0.21832
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.69318 0.8326
## Block (Intercept) 0.07665 0.2769
## Number of obs: 487, groups: obs, 487; Block, 3
##
## Fixed effects:
## Estimate Std. Error z value
## (Intercept) -2.01996 0.81642 -2.474
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 11.74179 1.30992 8.964
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -5.66290 0.68789 -8.232
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 1.03332 0.14306 7.223
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -0.06412 0.01022 -6.274
## Genetic_DiversityMedium -0.55984 0.10012 -5.592
## Genetic_DiversityLow -0.67366 0.09048 -7.445
## Pr(>|z|)
## (Intercept) 0.0134 *
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 < 2e-16 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 < 2e-16 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 5.08e-13 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 3.52e-10 ***
## Genetic_DiversityMedium 2.25e-08 ***
## Genetic_DiversityLow 9.69e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) s::(G_E,4,r=TRUE)1 s::(G_E,4,r=TRUE)2
## s::(G_E,4,r=TRUE)1 -0.965
## s::(G_E,4,r=TRUE)2 0.941 -0.993
## s::(G_E,4,r=TRUE)3 -0.915 0.976 -0.995
## s::(G_E,4,r=TRUE)4 0.889 -0.956 0.984
## Gntc_DvrstM -0.062 0.007 -0.008
## Gntc_DvrstL -0.065 0.004 -0.005
## s::(G_E,4,r=TRUE)3 s::(G_E,4,r=TRUE)4 Gnt_DM
## s::(G_E,4,r=TRUE)1
## s::(G_E,4,r=TRUE)2
## s::(G_E,4,r=TRUE)3
## s::(G_E,4,r=TRUE)4 -0.997
## Gntc_DvrstM 0.008 -0.009
## Gntc_DvrstL 0.005 -0.006 0.494
## convergence code: 0
## Model failed to converge with max|grad| = 0.053402 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 4.55 0.206 Inf 4.06 5.04
## Medium 3.93 0.216 Inf 3.42 4.45
## Low 3.69 0.201 Inf 3.21 4.17
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.619 0.187 Inf 3.310 0.0027
## High - Low 0.861 0.168 Inf 5.115 <.0001
## Medium - Low 0.242 0.184 Inf 1.319 0.3844
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates#######################
####################### MODEL INCLUDING ONLY RESCUED POPULATIONS
#######################
#Same but only for population not extinct
fit_fin <- glm(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity,
data = data[data$Extinction_finalstatus=="no",], family = "poisson")
filldata$Estimates_Withoutextpop <- predict(fit_fin, newdata=filldata, type = "response")
PLOT_Pop_size <- ggplot2::ggplot(data = data) +
geom_line(data = filldata, aes(x = Generation_Eggs,
y = Estimates_Withoutextpop,
colour = Genetic_Diversity), size = 1.5, linetype = "longdash") +
geom_point(data = data[data$Extinction_finalstatus=="no",],
aes(x = Generation_Eggs,
y = Nb_adults,
colour = Genetic_Diversity),
position = position_dodge(0.4), size = 1.4, alpha = 0.5) +
ylab("Population size") +
labs(color="Genetic diversity") +
xlab("Generation") +
theme_LO
PLOT_Pop_size
mod <- lme4::glmer(Nb_adults~stats::poly(Generation_Eggs,4,raw=TRUE)*Genetic_Diversity +
(1|Block) + (1|obs),
data = data[data$Extinction_finalstatus=="no",], family = "poisson")
mod1 <- lme4::glmer(Nb_adults~poly(Generation_Eggs,4,raw=TRUE)+Genetic_Diversity +
(1|Block) + (1|obs),
data = data[data$Extinction_finalstatus=="no",], family = "poisson")
anova(mod,mod1, test = "Chisq")## Data: data[data$Extinction_finalstatus == "no", ]
## Models:
## mod1: Nb_adults ~ poly(Generation_Eggs, 4, raw = TRUE) + Genetic_Diversity +
## mod1: (1 | Block) + (1 | obs)
## mod: Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +
## mod: (1 | Block) + (1 | obs)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 9 4759.9 4796.3 -2370.9 4741.9
## mod 17 4756.3 4825.1 -2361.1 4722.3 19.599 8 0.01196 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Nb_adults ~ stats::poly(Generation_Eggs, 4, raw = TRUE) * Genetic_Diversity +
## (1 | Block) + (1 | obs)
## Data: data[data$Extinction_finalstatus == "no", ]
##
## AIC BIC logLik deviance df.resid
## 4756.3 4825.1 -2361.1 4722.3 407
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.57063 -0.04879 0.00393 0.06967 0.29650
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs (Intercept) 0.36624 0.6052
## Block (Intercept) 0.02456 0.1567
## Number of obs: 424, groups: obs, 424; Block, 3
##
## Fixed effects:
## Estimate
## (Intercept) -1.9250178
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 10.7535583
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -5.0922588
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 0.9316535
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -0.0583188
## Genetic_DiversityMedium 1.1892488
## Genetic_DiversityLow 0.2794991
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium -2.0341925
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 1.0112012
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium -0.2039827
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.0137179
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow -0.4034287
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 0.0777534
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow -0.0101264
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.0006044
## Std. Error
## (Intercept) 0.9705855
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 1.5722049
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 0.8204596
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 0.1698215
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 0.0120893
## Genetic_DiversityMedium 1.5049889
## Genetic_DiversityLow 1.3816159
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 2.4448787
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 1.2751857
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.2640143
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.0188078
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 2.2495506
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 1.1751948
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.2435505
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.0173582
## z value
## (Intercept) -1.983
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 6.840
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 -6.207
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 5.486
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 -4.824
## Genetic_DiversityMedium 0.790
## Genetic_DiversityLow 0.202
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium -0.832
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 0.793
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium -0.773
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.729
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow -0.179
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 0.066
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow -0.042
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.035
## Pr(>|z|)
## (Intercept) 0.0473
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 7.93e-12
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 5.41e-10
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 4.11e-08
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 1.41e-06
## Genetic_DiversityMedium 0.4294
## Genetic_DiversityLow 0.8397
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium 0.4054
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium 0.4278
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium 0.4397
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium 0.4658
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow 0.8577
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow 0.9472
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow 0.9668
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow 0.9722
##
## (Intercept) *
## stats::poly(Generation_Eggs, 4, raw = TRUE)1 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)2 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)3 ***
## stats::poly(Generation_Eggs, 4, raw = TRUE)4 ***
## Genetic_DiversityMedium
## Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityMedium
## stats::poly(Generation_Eggs, 4, raw = TRUE)1:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)2:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)3:Genetic_DiversityLow
## stats::poly(Generation_Eggs, 4, raw = TRUE)4:Genetic_DiversityLow
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## convergence code: 0
## Model failed to converge with max|grad| = 3.25364 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
## - Rescale variables?
## Model is nearly unidentifiable: large eigenvalue ratio
## - Rescale variables?emmeans::emmeans(mod, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 4.53 0.133 Inf 4.21 4.85
## Medium 4.30 0.151 Inf 3.94 4.66
## Low 4.01 0.137 Inf 3.68 4.33
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.230 0.153 Inf 1.507 0.2874
## High - Low 0.523 0.140 Inf 3.746 0.0005
## Medium - Low 0.293 0.157 Inf 1.861 0.1503
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates#Signifcant interaction 3.4 Growth rate G6
head(data_phenotyping)## Population Week Block ID_Rep Divsersity Population.1 Box Initial.N Start
## 1 High1 5-week Block3 52 High 1 1 30 7/8/21
## 2 High1 5-week Block3 53 High 1 2 30 7/8/21
## 3 High13 5-week Block4 130 High 13 2 30 7/15/21
## 4 High13 5-week Block4 131 High 13 3 30 7/15/21
## 5 High13 5-week Block4 129 High 13 1 30 7/15/21
## 6 High14 5-week Block4 132 High 14 1 30 7/15/21
## End Larvae Pupae Adults Lambda Genetic_Diversity Nb_ttx p_adults
## 1 8/12/21 3 9 97 3.2333333 High 109 0.8899083
## 2 8/12/21 1 1 8 0.2666667 High 10 0.8000000
## 3 8/19/21 1 2 104 3.4666667 High 107 0.9719626
## 4 8/19/21 2 1 84 2.8000000 High 87 0.9655172
## 5 8/19/21 2 0 84 2.8000000 High 86 0.9767442
## 6 8/19/21 1 0 83 2.7666667 High 84 0.9880952
## p_pupae p_larvae He_remain He_remain_exp He_remain_end
## 1 0.08256881 0.027522936 0.9986008 0.9719424 0.9705825
## 2 0.10000000 0.100000000 0.9986008 0.9719424 0.9705825
## 3 0.01869159 0.009345794 0.9986008 0.9825008 0.9811261
## 4 0.01149425 0.022988506 0.9986008 0.9825008 0.9811261
## 5 0.00000000 0.023255814 0.9986008 0.9825008 0.9811261
## 6 0.00000000 0.011904762 0.9986008 0.9806907 0.9793186#Without considering extinct populations
#tapply(data_phenotyping$Lambda,data_phenotyping$Population,length)
#############################################################
################## Determine distribution ###################
#############################################################
## Poisson lognormal model
mlog <- lme4::lmer(log(Lambda+1) ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
#summary(mlog)
# Square root
msqrt <- lme4::lmer(sqrt(Lambda) ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
#summary(msqrt)
# Normal
mnorm <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
## Compare AIC
AIC(mlog, msqrt,mnorm)## df AIC
## mlog 7 169.8536
## msqrt 7 215.0823
## mnorm 7 626.7241#Square root a better fits to the data
#But we compare different values: Growth rate and Growth rate+1
# ## Simulate data
x.teo.log <- unlist(simulate(mlog))
x.teo.sqrt <- unlist(simulate(msqrt))
x.teo.norm <- unlist(simulate(mnorm))
# ## QQplot to compared Negative binomial and Poisson log normal distributions
qqplot(data_phenotyping$Lambda, x.teo.log,
main="QQ-plot", xlab="Observed data", ylab="Simulated data",
las=1, xlim = c(0, 10), ylim = c(0, 10), col='blue') ## QQ-plot
points(sort(data_phenotyping$Lambda), sort(x.teo.sqrt), pch=1, col='red')
points(sort(data_phenotyping$Lambda), sort(x.teo.norm), pch=1, col='green')
abline(0,1)
# ## Add legend
legend(0, 7, legend=c("Log", "Sqrt", "Normal"),
col=c("blue", "red", "green"),
pch=1, bty="n")
# ## Normal distribution provides a good fit to the data
# QQ plot
# qqnorm(resid(mlog))
# qqline(resid(mlog))
# qqnorm(resid(msqrt))
# qqline(resid(msqrt))
# qqnorm(resid(mnorm))
# qqline(resid(mnorm)) #best fit
#Distribution Residuals
hist(residuals(mlog),col="yellow",freq=F)
hist(residuals(msqrt),col="yellow",freq=F) #Seems the better fit
hist(residuals(mnorm),col="yellow",freq=F) #Seems the better fit
#Test normality
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest ## 1-mle-norm
## "not rejected"g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(msqrt), "norm"))
g1$kstest ## 1-mle-norm
## "not rejected"g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest## 1-mle-norm
## "not rejected"## Compare AIC
AIC(mlog, msqrt, mnorm)## df AIC
## mlog 7 169.8536
## msqrt 7 215.0823
## mnorm 7 626.7241#SLog a better fits to the data
#############################################################
################## Analysis ###################
#############################################################
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
summary(mod0)## Linear mixed model fit by REML ['lmerMod']
## Formula: Lambda ~ Genetic_Diversity + Block + (1 | Population)
## Data: data_phenotyping
##
## REML criterion at convergence: 612.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4322 -0.6449 -0.0938 0.5243 3.3448
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.2958 0.5439
## Residual 0.7591 0.8713
## Number of obs: 217, groups: Population, 77
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.5956 0.1904 13.629
## Genetic_DiversityMedium -1.1229 0.2271 -4.944
## Genetic_DiversityLow -0.8852 0.2128 -4.159
## BlockBlock4 0.2452 0.2059 1.191
## BlockBlock5 -0.1742 0.2407 -0.724
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM Gnt_DL BlckB4
## Gntc_DvrstM -0.488
## Gntc_DvrstL -0.559 0.401
## BlockBlock4 -0.605 0.044 0.075
## BlockBlock5 -0.593 0.111 0.185 0.457mod1 <- lme4::lmer(Lambda ~ Block + (1|Population),
data=data_phenotyping)
anova(mod0,mod1,test="Chisq") # difference## Data: data_phenotyping
## Models:
## mod1: Lambda ~ Block + (1 | Population)
## mod0: Lambda ~ Genetic_Diversity + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod1 5 641.07 657.97 -315.54 631.07
## mod0 7 618.35 642.01 -302.18 604.35 26.718 2 1.579e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
###### Posthoc
############
emmeans::emmeans(mod0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 2.62 0.136 52.2 2.28 2.96
## Medium 1.50 0.183 71.3 1.05 1.94
## Low 1.73 0.165 83.1 1.33 2.14
##
## Results are averaged over the levels of: Block
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 1.123 0.228 63.5 4.931 <.0001
## High - Low 0.885 0.213 69.5 4.147 0.0003
## Medium - Low -0.238 0.242 74.8 -0.983 0.5901
##
## Results are averaged over the levels of: Block
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates############
###### Predict
############
data_predict_extinct <- data.frame(Genetic_Diversity = levels(data_phenotyping$Genetic_Diversity),
Block = "Block4")
data_predict_extinct$predict <- predict(mod0,
type = "response",
re.form = NA,
newdata = data_predict_extinct)3.5 Remaining He difference
Rmisc::summarySE(data_G2,
measurevar = c("He_remain"),
groupvars = c("Genetic_Diversity"),
na.rm=TRUE)## Genetic_Diversity N He_remain sd se ci
## 1 High 27 0.9986008 0.000000000 0.000000000 0.000000000
## 2 Medium 23 0.6791246 0.006059307 0.001263453 0.002620241
## 3 Low 34 0.5441887 0.012687849 0.002175948 0.004427000mfull <- lm(He_remain ~ Genetic_Diversity, data=data_G2)
mless <- lm(He_remain ~ 1, data=data_G2)
anova(mfull,mless) #Keep effect## Analysis of Variance Table
##
## Model 1: He_remain ~ Genetic_Diversity
## Model 2: He_remain ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 81 0.0061
## 2 83 3.1868 -2 -3.1807 21048 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
###### Posthoc
############
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 0.9986 0.001673 81 0.9945 1.0027
## Medium 0.6791 0.001812 81 0.6747 0.6835
## Low 0.5442 0.001491 81 0.5406 0.5478
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.319 0.00247 81 129.527 <.0001
## High - Low 0.454 0.00224 81 202.800 <.0001
## Medium - Low 0.135 0.00235 81 57.498 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimatesRmisc::summarySE(data_G2,
measurevar = c("He_remain"),
groupvars = c("Genetic_Diversity"),
na.rm=TRUE)## Genetic_Diversity N He_remain sd se ci
## 1 High 27 0.9986008 0.000000000 0.000000000 0.000000000
## 2 Medium 23 0.6791246 0.006059307 0.001263453 0.002620241
## 3 Low 34 0.5441887 0.012687849 0.002175948 0.004427000#Calcul for end:
Rmisc::summarySE(data_G2[data_G2$Extinction_finalstatus!="yes",],
measurevar = c("He_remain_end"),
groupvars = c("Genetic_Diversity"),
na.rm=TRUE)## Genetic_Diversity N He_remain_end sd se ci
## 1 High 27 0.9697741 0.01868272 0.003595491 0.007390637
## 2 Medium 18 0.6482803 0.02453741 0.005783522 0.012202164
## 3 Low 26 0.5146307 0.02759467 0.005411760 0.011145728mfull <- lm(He_remain_end ~ Genetic_Diversity, data=data_G2[data_G2$Extinction_finalstatus!="yes",])
mless <- lm(He_remain_end ~ 1, data=data_G2[data_G2$Extinction_finalstatus!="yes",])
anova(mfull,mless) #Keep effect## Analysis of Variance Table
##
## Model 1: He_remain_end ~ Genetic_Diversity
## Model 2: He_remain_end ~ 1
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 68 0.03835
## 2 70 2.91180 -2 -2.8735 2547.7 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
###### Posthoc
############
emmeans::emmeans(mfull, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 0.970 0.00457 68 0.959 0.981
## Medium 0.648 0.00560 68 0.635 0.662
## Low 0.515 0.00466 68 0.503 0.526
##
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.321 0.00723 68 44.491 <.0001
## High - Low 0.455 0.00653 68 69.754 <.0001
## Medium - Low 0.134 0.00728 68 18.355 <.0001
##
## P value adjustment: tukey method for comparing a family of 3 estimates3.6 Growth rate and He
#Best model with genetic diversity
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping)
mod0 <- lme4::lmer(Lambda ~ Genetic_Diversity + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
mod1 <- lme4::lmer(Lambda ~ He_remain_end + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
mod2 <- lme4::lmer(Lambda ~ Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
mod3 <- lme4::lmer(Lambda ~ log(He_remain_end) + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
mod5 <- lme4::lmer(Lambda ~ exp(He_remain_end) + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
MuMIn::model.sel(mod0, mod1,mod2,mod3,mod4,mod5)## Model selection table
## (Int) Blc Gnt_Dvr He_rmn_end log(He_rmn_end) exp(He_rmn_end)
## mod3 2.6620 + 1.77
## mod1 0.2305 + 2.471
## mod5 -0.4063 + 1.151
## mod0 2.6630 + +
## mod2 1.9410 +
## mod4 1.9410 +
## family df logLik AICc delta weight
## mod3 gaussian(identity) 6 -283.053 578.5 0.00 0.473
## mod1 gaussian(identity) 6 -283.276 579.0 0.45 0.379
## mod5 gaussian(identity) 6 -284.363 581.2 2.62 0.128
## mod0 gaussian(identity) 7 -285.136 584.8 6.31 0.020
## mod2 gaussian(identity) 5 -297.218 604.7 26.21 0.000
## mod4 gaussian(identity) 5 -297.218 604.7 26.21 0.000
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | Population'#Makes sense: best model log He
x <- seq(1:100)
y <- seq(1001:1100)
plot(log(x),y)
PLOT_He_withextinct
#############################################################
################## Analysis ###################
#############################################################
mod3 <- lme4::lmer(Lambda ~ log(He_remain_end) + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
mod4 <- lme4::lmer(Lambda ~ 1 + Block + (1|Population),
data=data_phenotyping[!is.na(data_phenotyping$He_remain_end),])
anova(mod3,mod4,test="Chisq") # difference## Data: data_phenotyping[!is.na(data_phenotyping$He_remain_end), ]
## Models:
## mod4: Lambda ~ 1 + Block + (1 | Population)
## mod3: Lambda ~ log(He_remain_end) + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## mod4 5 599.98 616.53 -294.99 589.98
## mod3 6 571.83 591.68 -279.91 559.83 30.158 1 3.983e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(mod3)## Linear mixed model fit by REML ['lmerMod']
## Formula: Lambda ~ log(He_remain_end) + Block + (1 | Population)
## Data: data_phenotyping[!is.na(data_phenotyping$He_remain_end), ]
##
## REML criterion at convergence: 566.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4098 -0.6488 -0.0811 0.5245 3.3605
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.2616 0.5115
## Residual 0.7609 0.8723
## Number of obs: 202, groups: Population, 73
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 2.6623 0.1900 14.015
## log(He_remain_end) 1.7700 0.3028 5.846
## BlockBlock4 0.2491 0.2060 1.209
## BlockBlock5 -0.1118 0.2391 -0.468
##
## Correlation of Fixed Effects:
## (Intr) l(H__) BlckB4
## lg(H_rmn_n) 0.634
## BlockBlock4 -0.617 -0.103
## BlockBlock5 -0.569 -0.148 0.454#############################################################
################## Predict ###################
#############################################################
data_predict_lambda <- data.frame(He_remain_end =
seq(min(data_phenotyping[!is.na(data_phenotyping$He_remain_end),]$He_remain_end),
max(data_phenotyping[!is.na(data_phenotyping$He_remain_end),]$He_remain_end),
0.01),
Block = "Block4")
data_predict_lambda$lambda_predict <- predict(mod3,
type = "response",
re.form = NA,
newdata = data_predict_lambda)
PLOT_He_expect <- ggplot2::ggplot() +
geom_line(data = data_predict_lambda,
aes(x = He_remain_end, y = lambda_predict),
linetype = "longdash", col = "grey40", size = 1.25) +
geom_point(data = SUM_POP_He_ALL_extinct,
aes(x = He_remain_end, y = Lambda,
colour = Genetic_Diversity,
shape = Extinction_finalstatus,
size = N), alpha = 0.7) +
ylab("Growth rate") +
xlab("Remaining heterozygosity") +
scale_shape_manual(values = c(19, 1)) +
theme_LO +
labs(color="Genetic diversity") +
labs(shape="Extinct population") +
labs(size="Replicates")
PLOT_He_expect
#
# cowplot::save_plot(here::here("figures","Fitness_He_final.pdf"), PLOT_He_withextinct,
# base_height = 12/cm(1), base_width = 18/cm(1), dpi = 610)3.7 Adults G6
m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Population) + (1|ID_Rep:Population),
family = "poisson", data = data_phenotyping)
summary(m0)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
## Data: data_phenotyping
##
## AIC BIC logLik deviance df.resid
## 2154.4 2171.3 -1072.2 2144.4 212
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1441 -0.1230 0.0254 0.1360 0.3923
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_Rep:Population (Intercept) 0.2041 0.4518
## Population (Intercept) 0.2807 0.5298
## Number of obs: 217, groups: ID_Rep:Population, 217; Population, 77
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.2792 0.1105 38.715 < 2e-16 ***
## Genetic_DiversityMedium -0.7386 0.1812 -4.076 4.59e-05 ***
## Genetic_DiversityLow -0.5250 0.1664 -3.156 0.0016 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM
## Gntc_DvrstM -0.611
## Gntc_DvrstL -0.667 0.413# Equivalent to:
# m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Population/ID_Rep),
# family = "poisson", data = data_5week)
#dispersion_lme4::glmer(m0) #dispersion ok
#Note: (1|School) + (1|Class:School) is equivalent to (1|School/Class)
m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Population) + (1|ID_Rep:Population),
family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1 + (1|Population) + (1|ID_Rep:Population),
family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") # difference## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Population) + (1 | ID_Rep:Population)
## m0: Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 3 2167.5 2177.6 -1080.7 2161.5
## m0 5 2154.4 2171.3 -1072.2 2144.4 17.04 2 0.0001994 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1summary(m0)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula:
## Adults ~ Genetic_Diversity + (1 | Population) + (1 | ID_Rep:Population)
## Data: data_phenotyping
##
## AIC BIC logLik deviance df.resid
## 2154.4 2171.3 -1072.2 2144.4 212
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1441 -0.1230 0.0254 0.1360 0.3923
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID_Rep:Population (Intercept) 0.2041 0.4518
## Population (Intercept) 0.2807 0.5298
## Number of obs: 217, groups: ID_Rep:Population, 217; Population, 77
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.2792 0.1105 38.715 < 2e-16 ***
## Genetic_DiversityMedium -0.7386 0.1812 -4.076 4.59e-05 ***
## Genetic_DiversityLow -0.5250 0.1664 -3.156 0.0016 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM
## Gntc_DvrstM -0.611
## Gntc_DvrstL -0.667 0.413############
###### Posthoc
############
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 4.28 0.111 Inf 4.02 4.54
## Medium 3.54 0.143 Inf 3.20 3.88
## Low 3.75 0.124 Inf 3.46 4.05
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.739 0.181 Inf 4.076 0.0001
## High - Low 0.525 0.166 Inf 3.156 0.0046
## Medium - Low -0.214 0.189 Inf -1.132 0.4943
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates#test double nested
m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Genetic_Diversity/Population/ID_Rep),
family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1 + (1|Genetic_Diversity/Population/ID_Rep),
family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") ## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/Population/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/Population/ID_Rep)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 4 2160.9 2174.4 -1076.4 2152.9
## m0 6 2156.4 2176.7 -1072.2 2144.4 8.4441 2 0.01467 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1m0 <- lme4::glmer(Adults ~ Genetic_Diversity + (1|Genetic_Diversity/ID_Rep),
family = "poisson", data = data_phenotyping)
m1 <- lme4::glmer(Adults ~ 1 + (1|Genetic_Diversity/ID_Rep),
family = "poisson", data = data_phenotyping)
anova(m0,m1,test="Chisq") ## Data: data_phenotyping
## Models:
## m1: Adults ~ 1 + (1 | Genetic_Diversity/ID_Rep)
## m0: Adults ~ Genetic_Diversity + (1 | Genetic_Diversity/ID_Rep)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 3 2189.9 2200.0 -1091.9 2183.9
## m0 5 2183.0 2199.9 -1086.5 2173.0 10.858 2 0.004387 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 13.8 Adults G2
m0 <- glm(Nb_adults ~ Genetic_Diversity + Block, family = "poisson", data = data_G2)
summary(m0)##
## Call:
## glm(formula = Nb_adults ~ Genetic_Diversity + Block, family = "poisson",
## data = data_G2)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -18.7057 -2.7621 0.2126 3.3983 11.4454
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.73890 0.01497 383.312 < 2e-16 ***
## Genetic_DiversityMedium -0.19408 0.01628 -11.920 < 2e-16 ***
## Genetic_DiversityLow -0.19278 0.01460 -13.205 < 2e-16 ***
## BlockBlock4 0.10767 0.01579 6.817 9.3e-12 ***
## BlockBlock5 0.17743 0.01600 11.088 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 2379.6 on 83 degrees of freedom
## Residual deviance: 2027.9 on 79 degrees of freedom
## AIC: 2667.2
##
## Number of Fisher Scoring iterations: 4m1<- glm(Nb_adults ~ Block, family = "poisson", data = data_G2)
anova(m0,m1,test="Chisq") # difference## Analysis of Deviance Table
##
## Model 1: Nb_adults ~ Genetic_Diversity + Block
## Model 2: Nb_adults ~ Block
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 79 2027.9
## 2 81 2241.4 -2 -213.52 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
###### Posthoc
############
emmeans::emmeans(m0, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df asymp.LCL asymp.UCL
## High 5.834 0.01051 Inf 5.809 5.859
## Medium 5.640 0.01243 Inf 5.610 5.670
## Low 5.641 0.01022 Inf 5.617 5.666
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df z.ratio p.value
## High - Medium 0.1941 0.0163 Inf 11.920 <.0001
## High - Low 0.1928 0.0146 Inf 13.205 <.0001
## Medium - Low -0.0013 0.0161 Inf -0.081 0.9964
##
## Results are averaged over the levels of: Block
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates3.9 Evolution of growthrate over time
############
###### Distribution
############
hist(data_evolution_Lambda$Lambda, breaks=50)
# Square root
msqrt <- lme4::lmer(sqrt(Lambda) ~ Genetic_Diversity*Phenotyping +
Block + (1|Population),
data = data_evolution_Lambda)
summary(msqrt)## Linear mixed model fit by REML ['lmerMod']
## Formula: sqrt(Lambda) ~ Genetic_Diversity * Phenotyping + Block + (1 |
## Population)
## Data: data_evolution_Lambda
##
## REML criterion at convergence: 252.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4979 -0.4737 0.0672 0.5620 2.4161
##
## Random effects:
## Groups Name Variance Std.Dev.
## Population (Intercept) 0.0295 0.1718
## Residual 0.1040 0.3225
## Number of obs: 301, groups: Population, 88
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 1.779696 0.081683 21.788
## Genetic_DiversityMedium -0.171161 0.103522 -1.653
## Genetic_DiversityLow -0.150003 0.094035 -1.595
## PhenotypingFinal -0.243939 0.070714 -3.450
## BlockBlock4 0.102764 0.063883 1.609
## BlockBlock5 0.003821 0.070206 0.054
## Genetic_DiversityMedium:PhenotypingFinal -0.251050 0.110515 -2.272
## Genetic_DiversityLow:PhenotypingFinal -0.154999 0.101256 -1.531
##
## Correlation of Fixed Effects:
## (Intr) Gnt_DM Gnt_DL PhntyF BlckB4 BlckB5 G_DM:P
## Gntc_DvrstM -0.600
## Gntc_DvrstL -0.655 0.506
## PhntypngFnl -0.685 0.534 0.588
## BlockBlock4 -0.462 0.054 0.041 0.028
## BlockBlock5 -0.424 0.006 0.009 0.002 0.488
## Gntc_DvM:PF 0.428 -0.745 -0.377 -0.641 -0.028 0.050
## Gntc_DvL:PF 0.444 -0.373 -0.733 -0.699 0.012 0.088 0.453# log
mlog <- lme4::lmer(log(Lambda+1) ~ Genetic_Diversity*Phenotyping +
Block + (1|Population),
data=data_evolution_Lambda)
# Normal
mnorm <- lme4::lmer(Lambda ~ Genetic_Diversity*Phenotyping +
Block + (1|Population),
data=data_evolution_Lambda)
## Compare AIC
AIC(mlog, msqrt,mnorm)## df AIC
## mlog 10 207.0029
## msqrt 10 272.0671
## mnorm 10 848.9685#Square root a better fits to the data
#But we compare different values: Growth rate and Growth rate+1
#Distribution Residuals
hist(residuals(mlog),col="yellow",freq=F)
hist(residuals(msqrt),col="yellow",freq=F) #Seems the better fit
hist(residuals(mnorm),col="yellow",freq=F)
#Test normality
g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mlog), "norm"))
g1$kstest ## 1-mle-norm
## "rejected"g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(msqrt), "norm"))
g1$kstest ## 1-mle-norm
## "not rejected"g1 <- fitdistrplus::gofstat(fitdistrplus::fitdist(residuals(mnorm), "norm"))
g1$kstest## 1-mle-norm
## "not rejected"# QQ plot
# plot(mlog, which = 2)
# plot(msqrt, which = 2) #Seems the better fit!
# plot(mnorm, which = 2) #Seems the better fit but very high AIC
qqnorm(resid(mlog))
qqline(resid(mlog))
qqnorm(resid(msqrt))
qqline(resid(msqrt))
qqnorm(resid(mnorm))
qqline(resid(mnorm)) #best fit
############
###### Model selection
############
m0 <- lme4::lmer(Lambda ~ Genetic_Diversity*Phenotyping +
Block + (1|Population) ,
data=data_evolution_Lambda)
m1 <- lme4::lmer(Lambda ~ Genetic_Diversity+Phenotyping +
Block + (1|Population) ,
data=data_evolution_Lambda)
anova(m0,m1,test="Chisq") #Remove interaction## Data: data_evolution_Lambda
## Models:
## m1: Lambda ~ Genetic_Diversity + Phenotyping + Block + (1 | Population)
## m0: Lambda ~ Genetic_Diversity * Phenotyping + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m1 8 832.93 862.59 -408.46 816.93
## m0 10 834.25 871.32 -407.12 814.25 2.6798 2 0.2619#Same result with:
# m1 <- lmer(Lambda ~ Genetic_Diversity+Phenotyping +
# (1|Block) + (1|Population) + (1|ID_Rep:Population),
# data=data_evolution_Lambda)
# OR
# m1 <- lmer(Lambda ~ Genetic_Diversity+Phenotyping +
# (1|Block) + (1|Population) ,
# data=data_evolution_Lambda)
m2 <- lme4::lmer(Lambda ~ Phenotyping +
Block + (1|Population),
data=data_evolution_Lambda)
m3 <- lme4::lmer(Lambda ~ Genetic_Diversity +
Block + (1|Population),
data=data_evolution_Lambda)
anova(m1,m2,test="Chisq") #Keep phenotyping## Data: data_evolution_Lambda
## Models:
## m2: Lambda ~ Phenotyping + Block + (1 | Population)
## m1: Lambda ~ Genetic_Diversity + Phenotyping + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m2 6 855.30 877.54 -421.65 843.30
## m1 8 832.93 862.59 -408.46 816.93 26.369 2 1.88e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m1,m3,test="Chisq") #Keep genetic diversity## Data: data_evolution_Lambda
## Models:
## m3: Lambda ~ Genetic_Diversity + Block + (1 | Population)
## m1: Lambda ~ Genetic_Diversity + Phenotyping + Block + (1 | Population)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## m3 7 896.49 922.44 -441.24 882.49
## m1 8 832.93 862.59 -408.46 816.93 65.558 1 5.642e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1############
###### Posthoc
############
mod_final <- lme4::lmer(Lambda ~ Genetic_Diversity+Phenotyping +
Block + (1|Population),
data=data_evolution_Lambda)
emmeans::emmeans(mod_final, list(pairwise ~ Genetic_Diversity), adjust = "tukey") ## $`emmeans of Genetic_Diversity`
## Genetic_Diversity emmean SE df lower.CL upper.CL
## High 3.05 0.120 69.4 2.76 3.34
## Medium 2.14 0.143 87.6 1.79 2.48
## Low 2.32 0.123 99.8 2.02 2.62
##
## Results are averaged over the levels of: Phenotyping, Block
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 3 estimates
##
## $`pairwise differences of Genetic_Diversity`
## contrast estimate SE df t.ratio p.value
## High - Medium 0.916 0.184 72.8 4.988 <.0001
## High - Low 0.734 0.170 78.2 4.316 0.0001
## Medium - Low -0.182 0.187 87.5 -0.972 0.5963
##
## Results are averaged over the levels of: Phenotyping, Block
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimatesemmeans::emmeans(mod_final, list(pairwise ~ Phenotyping), adjust = "tukey") ### $`emmeans of Phenotyping`
## Phenotyping emmean SE df lower.CL upper.CL
## Initial 3.00 0.1072 254 2.76 3.24
## Final 2.01 0.0837 109 1.82 2.20
##
## Results are averaged over the levels of: Genetic_Diversity, Block
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## Conf-level adjustment: sidak method for 2 estimates
##
## $`pairwise differences of Phenotyping`
## contrast estimate SE df t.ratio p.value
## Initial - Final 0.993 0.117 250 8.474 <.0001
##
## Results are averaged over the levels of: Genetic_Diversity, Block
## Degrees-of-freedom method: kenward-roger#